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Preface 


This  thesis  describes  a  method  used  at  the  Air  Force  Test  Pilot 
School  of  identifying  aircraft  parameters  via  frequency  response 
analysis  techniques.  This  method  uses  various  input  and  response 
signals  from  the  aircraft  and  transforms  these  signals  to  generate 
frequency  response  data  (Bode  Plots).  Lower  order  equivalent  systems 
were  fit  to  the  frequency  response  data  and  equivalent  aircraft 
parameters  were  determined.  The  method  was  developed  using  simulated 
aircraft  time  histories  which  were  analytically  derived.  Flight 
testing  was  later  accomplished  which  verified  and  supplemented  the 
analysis . 

The  work  on  this  thesis  was  accomplished  through  the  Joint  Air 
Force  Institute  of  Technology /Test  Pilot  School  (AFIT/TPS)  Program. 

I  wish  to  express  my  gratitude  to  my  thesis  advisors  Major  (Dr.) 
James  T.  Silverthorn  and  Dr.  Robert  A.  Calico.  Much  credit  also  is 
due  the  USAF  Test  Pilot  School  for  sponsoring  a  test  program  through 
which  my  analytic  work  could  be  tested  and  verified.  I  also  wish  to 
thank  my  wife,  Deborah,  for  her  support  throughout  the  two  year  joint 
AFIT/TPS  program. 


Allan  T.  Reed 
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Abst  ract 


This  paper  presents  the  results  of  a  project  which  used  frequency- 
analysis  methods  applied  to  flight  test  data  in  order  to  identify 
aircraft  parameters.  Computer  programs  were  developed  to  generate 
simulated  flight  test  data  so  the  frequency  response  programs  could  be 
tested  using  a  noise  free  data  source.  Once  the  simulated  data 
programs  were  complete,  the  frequency  response  methods  were  developed. 
The  frequency  response  method  uses  the  cross-spectral  density 
technique  to  generate  magnitude  and  phase  information.  The  program 
also  generates  the  power  spectral  density  functions.  Noise 
contamination  studies  were  made  and  the  frequency  response  program  was 
modified  to  use  a  window-averaging  technique  to  reduce  the  effects  of 
noise  as  well  as  to  generate  a  coherence  function  to  display  where 
poor  correlation  between  input  and  output  was  occurring  due  to  noise. 
The  programs  were  used  with  the  test  project,  HAVE  DELAY,  which  was 
conducted  at  the  USAF  Test  Pilot  School.  Lower  order  equivalent 
systems  techniques  were  used  to  fit  results  from  flight  test  to  lower 
order  systems.  The  frequency  response  program  worked  well  up  to  a 
frequency  of  about  15  radians  per  second.  Above  15  radians  per  second 
the  program  suffers  from  additional  noise  and  aliasing  effects. 
Recommendations  were  made  to  study  means  of  reducing  noise  effects  to 
include  anti-aliasing  filters  and  noise  processing  schemes  for  digital 
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EVALUATION  OF  A  FREQUENCY  RESPONSE  TECHNIQUE 


FOR  AIRCRAFT  SYSTEM  IDENTIFICATION 


I.  Introduction 


This  thesis  describes  a  method  of  identifying  aircraft  dynamic 
parameters  (natural  frequencies,  damping  ratios,  etc.)  using  a 
frequency  response  technique.  These  techniques  were  developed  for  use 
by  the  USAF  Test  Pilot  School. 

Background 

As  modern  aircraft  and  the  flight  control  systems  become  more 
advanced,  the  response  of  these  aircraft  become  increasingly  difficult 
to  evaluate.  Most  modern  aircraft  no  longer  exhibit  the  classic 
response  characteristics  for  which  exists  a  large  body  of  knowledge 
(i.e.  short  period,  phugoid,  dutch  roll,  roll,  and  spiral).  This  has 
made  the  determination  of  MIL-SPEC  (Ref  3),  compliance  extremely 
difficult,  if  not  impossible.  Also  many  times  the  design  engineer  is 
more  interested  in  how  the  aircraft  responds  during  pilot-in-the-loop 
or  operational  type  maneuvering  than  how  that  aircraft  responds  open 
loop. 

When  flight  control  systems  played  a  minor  role  in  determining 
overall  aircraft  response  then  open  loop  response  was  adequate.  But 
with  complex  flight  control  systems  which  significantly  alter  the 
nature  of  the  response,  the  flight  test  engineer  needs  to  examine  the 


total  system  response— that  is,  from  the  force  the  pilot  applies  to 


the  control  stick  to  the  final  aircraft  response.  Often  only  part  of 
the  system  needs  to  be  analyzed. 

A  frequency  response  analyzer  is  ideally  suited  for  this  problem. 
Unlike  classical  flight  test  techniques  which  require  precise  inputs 
and  test  parameters  to  be  performed  by  the  test  pilot,  frequency 
analysis  methods  only  requires  that  the  system  be  sufficiently  excited 
for  a  sufficient  amount  of  time  (about  20-40  sec),  over  a  sufficient 
bandwidth  (about  10  rad/sec). 

Frequency  response  analysis  methods  are  also  superior  to 
classical  techniques  in  that  frequency  response  techniques  can  be  used 
while  the  pilot  is  flying  tasks  which  are  not  much  different  from  tasks 
the  pilot  would  be  flying  if  he  were  employing  the  aircraft 
operationally.  Typical  tasks  for  frequency  analysis  techniques  would 
be  air  refueling,  air-to-air  tracking,  air-to-ground  tracking, 
precision  formation,  or  precision  landing.  Thus,  the  design  engineer 
is  analyzing  the  total  aircraft  system  (airframe  plus  flight 
controls),  and  is  doing  so  with  the  aircraft  flying  maneuvers  similar 
to  those  flown  operationally. 

Current  Flight  Test  Techniques 

The  United  States  Air  Force  Test  Pilot  School  currently  uses  open 
loop  testing  to  determine  aircraft  natural  frequencies  and  damping 
ratios.  The  existing  technique  involve  stabilizing  the  aircraft  at 
the  desired  test  conditions  (trim  shot),  and  then  applying  a  test 
input,  (release  from  sideslip,  doublet,  etc.),  releasing  controls,  and 
then  recording  aircraft  response.  Obviously  the  response  being 
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evaluated  is  an  open  loop  response  and  effects  of  the  flight  controls 


are  not  entirely  present.  Data  reduction  for  this  technique  involves 
the  examination  of  strip  charts  (if  the  aircraft  is  instrumented)  to 
determine  the  natural  frequency  and  damping  ratios.  Damping  ratios 
are  typically  determined  by  a  time  ratio  method  or  log-decrement 
technique.  (Noninstrumented  aircraft  are  evaluated  by  the  pilot 
timing  and  counting  the  number  of  overshoots  following  the  input  to 
determine  natural  frequency  and  damping  ratio). 

There  are  several  problems  associated  with  the  current 
techniques.  Aircraft  response  is  generally  a  combination  of  motions. 
In  other  words,  motion  about  the  pitch  axis  is  virtually  never  just 
short  period  motion  or  just  phugoid  motion,  but  rather  some 
combination  of  the  two  motions.  Only  by  exciting  an  aircraft  by 
introducing  as  initial  conditions  a  motion  which  is  some  multiple  of 
the  desired  eiqenvector  can  one  get  motion  of  a  single  mode  as  a 
response.  Often  a  sinusoidal  pulse  or  control  doublet  is  performed  by 
the  pilot  at  a  frequency  which  is  very  close  to  the  natural  frequency 
of  the  mode  being  analyzed.  The  current  methods  and  frequency  sweeps 
can  provide  an  idea  of  the  frequency  response  based  upon  flight  test 
and  data  reduction.  Classical  techniques  tell  much  about  an  aircraft. 
But,  they  are  free  response  techniques — they  tell  what  happens  to  the 
aircraft  when  disturbed  from  equilibrium,  controls  released,  and 
resulting  motion  taking  place.  What  really  needs  to  be  examined  is 
the  forced  response — because  the  critical  phases  of  flight  such  as 
air-to-air,  air-to-ground,  air  refueling,  formation,  takeoff,  and 
landing  are  forced  responses  (pilot-in-the-loop).  Furthermore,  modern 
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aircraft  are  often  so  augmented  with  feedback  for  stability  and 
control  that  the  classic  modes  may  not  appear  in  the  standard  manner 
or  may  be  surrounded  by  additional  modes  introduced  by  flight  controls 
(feel  and  trim  systems,  etc.)*  As  an  example  of  inherent  errors  in 
the  classic  methods  of  exciting  an  aircraft  in-flight,  consider  the 
case  where  the  engineer  wishes  to  record  the  free  aircraft  response  to 
a  doublet  input  and  typical  pilot  realization  of  that  elevator  input 
as  shown  in  Figure  1.  Shown  are  the  elevator  input  signals  as 
typically  flown.  Obviously  a  discrepancy  exists  between  the  desired 
input  and  the  resulting  input,  and  therefore,  uncertainties  will  exist 
as  to  the  real  significance  of  the  recorded  response  to  these  signals. 

The  concept  of  examining  the  forced  response  of  an  aircraft  is 
important  as  forced  response  is  the  basis  for  handling  qualities 
evaluations.  The  forced  response,  also  called  mission  oriented 
response  or  pilot-in-the-loop  response  (Ref  14),  is  what  ultimately 
must  be  optimized  in  the  design  of  flight  control  systems.  Handling 
qualities  discovered  during  evaluation  of  this  forced  response  are  a 
function  of  the  total  system,  from  the  pilot  applying  stick  inputs 
down  to  the  airframe  aerodynamics. 

The  aeronautical/control  engineer  seeks  better  methods  to  analyze 
and  optimize  the  forced  response  of  an  aircraft.  This  thesis  presents 
a  method  of  analyzing  aircraft  via  frequency  response  methods. 
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DESIRED  DOUBLET 


DESIRED  ELEVATOR  DEFLECTION  HISTORY 


ELEVATOR 

DEFLECTION 


ELEVATOR 

DEFLECTION 


Figure  1  Desired  and  Typical  Elevator  Inputs 


II.  Frequency  Response  Method 


Frequency  response  methods  require  obtaining  time  histories  of 
various  inputs  and  the  respective  aircraft  response.  If  these  time 
histories  can  be  appropriately  transformed  into  the  frequency  domain 
then  classic  frequency  domain  analysis  methods  (i.e.,  Bode,  amplitude 
and  phase  plots)  serve  the  engineer  in  the  analysis  and  optimization 
process.  This  concept  is  shown  in  Figure  2. 

The  choice  of  transformations  is  logically  the  Fourier  Transform 
as  shown  in  Figure  3,  and  since  this  project  used  flight  test  data 
which  was  digitally  sampled,  the  transform  which  was  used  was  a 
Discrete  Fourier  Transform  or  DFT. 

Input  Signal 

The  choice  of  input  signals  is  perhaps  unclear.  If  one  demands  a 
rigorous,  precisely  flown  input  signal,  x(t),  then  the  same  problem  as 
the  open  loop  method  is  encountered.  The  problem  with  the  open  loop 
test  was  that  the  pilot  can  never  input  precisely  the  desired  signal 
(the  test  pilot  cannot  input  a  perfect  step,  a  perfect  doublet,  etc.) 
(Figure  1) 

Examination  of  the  concepts  behind  the  properties  of  the  Impulse 
function  reveal  an  interesting  property:  if  the  input  and  response 
time  signals  are  measured  simultaneously  and  then  Fourier  Transformed, 
the  resulting  division  of  these  two  signals,  i.e.,  Y(a>)/X(w) 
represents  the  transfer  function  between  input  and  output.  Thus, 
after  a  processing  of  test  data  from  a  single  sampling  of  input  and 
output  a  complete  frequency  response  plot  can  be  constructed.  The 
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t  s  time 

*»>  *  frequency 

h(t)  s  impulse  resp'nse  function 
H(w)  =  frequency  response  function 

Figure  2  Frequency  Analysis  Concept 


7 


WAVEFORM  DEFINED 


only  restriction  to  the  input  signal  is  that  the  frequency  content 
must  extend  across  the  desired  frequency  range  of  interest.  For 
aircraft  flight  dynamics,  control  engineers  are  most  concerned  with 
the  frequency  range  0.1  to  10  rad/sec. 

The  fact  that  the  input  signal  must  contain  sufficient  power  in 
the  frequency  range  of  interest  can  be  complied  without  any  additional 
restrictions  to  the  test  pilot's  input  signal — thus,  random  input 
signals  can  be  used  if  they  contain  an  adequate  amount  of  power  in  the 
range  0.1  to  10  rad/sec. 

Thus  the  approach  for  this  project  was  to  digitally  sample 
various  inputs  of  interest  and  to  simultaneous ly  sample  various 
outputs,  use  a  digital  Fourier  Transform  to  convert  signals  from  the 
time  domain  into  the  frequency  domain  and  then,  to  compute  and  plot  the 
frequency  response  transfer  function  magnitude  and  phase  angle  for 
analysis . 

Aliasing 

Several  problems  arise  in  the  frequency  analysis  of  sampled  data. 
One  such  problem  is  that  of  aliasing.  Consider  the  two  continuous 
waveforms  shown  in  Figure  A:  one  of  low  frequency  and  one  of  high 
frequency.  As  long  as  the  signal  is  sampled  at  least  as  fast  as  the 
theoretical  limit  (a  minimum  of  two  samples  per  cycle),  the  transform 
methods  will  correctly  identify  the  signals.  However,  there  are 
always  practical  limits  as  to  how  fast  the  signal  can  be  sampled 
(i.e.,  8  samples/sec,  16  samples/sec,  etc.).  A  problem  occurs  if 


there  are  frequencies  present  in  the  signals  which  are  higher  than  the 
maximum  frequency  which  can  be  identified  by  the  sampling  rate. 
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CONDITIONS  FOR  ALIASING 

nal  high  frequency  signal 

low  sample  rate  (large  sampling  interval) 
less  than  2  samples  per  cycle 


Figure  4  shows  an  example  of  the  classic  case  of  aliasing.  The  dotted 
line  shown  is  an  alias  of  the  signal  represented  by  the  solid  line. 

This  has  a  significant  impact  on  the  analysis  of  aircraft  data. 
Given  the  frequency  range  of  interest,  one  must  ensure  that  the 
sampling  rate  is  such  that  the  highest  frequency  is  sampled  at  least 
twice  per  cycle  (the  Nyquist  rate).  Aircraft  flight  dynamics 
generally  deal  with  a  frequency  range  of  0.1  to  10  rad/sec.  Therefore 
a  sampling  rate  is  needed  which  gives  at  least  two  samples  per  cycle 
of  a  10  rad/sec  signal,  i.e.  sample  at  least  twice  during  a  cycle  of 
0.63  seconds  time  interval  or  once  every  third  of  a  second.  Thus,  the 
bare  minimum  sample  rate  is  three  tmes  a  second  (Fn  *  2TT/2AT). 

This  is,  however,  a  theoretical  limit.  By  following  proven 
engineering  guidelines  of  five  samples  per  cycle,  one  finds  that  one 
should  sample  the  signals  at  approximately  0.13  sec  (Fn  *  2TT/5*2&X). 

Sixteen  samples  per  second  should  allow  good  analysis  out  to 
about  20  rad/sec.  Some  Test  Pilot  School  aircraft  have  data 
acquisition  systems  (DAS)  which  sample  20  times  per  second  which 
should  allow  analysis  out  to  about  25  rad/sec.  However,  remember  this 
discussion  of  aliasing  assumes  that  one  knows  the  highest  frequency 
present  in  the  data  and  can  sample  at  an  appropriate  rate  to  prevent 
the  higher  frequencies  from  aliasing  the  lower  frequencies.  Since  a 
basic  linear  system  is  assumed,  the  output  spectrum  will  be  limited  to 
the  same  range  as  the  input  spectrum. 

Once  the  proper  choice  of  sampling  rate  is  used,  the  aliasing 
will  be  reduced  to  effects  which  occur  from  noise  in  the  data.  If 


analysis  of  input  and  output  spectrum  show  aliasing  still  to  be  a 


problem  (i.e.,  due  Co  nonlinear  effects,  buffetting,  aerodynamic 
noise,  etc.)  then  a  low  pass  anti-aliasing  filter  and  higher  sample 
rate  should  reduce  aliasing  effects. 

Window  Effects 

Another  problem  which  arises  in  the  method  of  frequency 
transformations  of  sampled  data  is  the  fact  the  input  and  output 
signals  are  not  periodic  (we  are  dealing  with  "random"  inputs).  The 
fact  the  signals  are  not  periodic  cause  an  effect  known  as  "leakage". 

If  20  sec  of  data  are  analyzed,  the  effect  is  as  though  one  is 
looking  at  data  through  a  20  sec  long  rectangular  window  and 
neglecting  everything  which  happened  before  the  window  started  and 
everything  which  happened  after  the  20  sec  window  stopped. 

In  essence,  multiplication  of  the  data  with  the  window  in  the 
time  domain  produce  a  convolution  of  the  Fourier  transform  of  the  data 
window  W(u),  with  the  Fourier  transform  of  the  data  X(*>). 
Unfortunately,  the  rectangular  data  window  has  many  undesirable 
features  in  the  frequency  domain.  Namely,  its  flat  top  shape  results 
in  a  multilobe  shape  in  the  frequency  domain.  This  multilobe  shape 
will  allow  energy  to  appear  at  these  lobes  (side  lobe),  thereby  giving 
the  impression  that  nearby  frequencies  are  present  when  in  fact  they 
are  not. 

Consider  the  examples  shown  in  Figures  5,  6,  and  7.  Figures  5 
and  6  show  periodic  and  non-peroidic  signals.  Figure  7  shows  a  signal 
which  is  periodic  in  the  window  0-T  sec.  Upon  transformation  to  the 
frequency  domain,  the  positive  frequencies  would  appear  as  shown. 

As  expected,  the  transformation  of  the  2sin3t  waveform  into  the 


frequency  domain  shows  chat  the  only  signal  present  is  at  3  rad/sec. 
Now  consider  the  signal  of  Figure  8.  This  signal  is  not  periodic  in 
the  T  sec  window.  Notice  that  the  3  rad/sec  signal  does  not  appear  as 
a  single  spike  at  3  rad/sec  but  rather  is  smeared  from  about  2.5  to 
3.5  rad/sec.  Thus  "leakage"  has  occurred  through  the  side  lobe  of  the 
data  window  w(t)  because  the  data  signal  was  not  periodic  in  the 
window. 

In  aircraft  frequency  response  analysis,  virtually  all  of  the 
signals  (input  and  output)  will  be  non-periodic,  and  therefore,  one 
must  account  for  leakage  or  else  the  frequency  spectrum  will  be  so 
smeared  as  to  make  frequency  analysis  impossible.  The  only  other  cure 
for  the  leakage  problem  is  to  take  an  extremely  long  sample  size — with 
the  limiting  case  being  no  leakage  with  an  infinitely  long  window  or 
sample  size. 

Window  Selection 

The  practical  way  to  cure  the  leakage  problem  is  in  the  selection 
of  a  data  window  other  than  a  rectangular  window.  The  method  of 
choosing  the  best  window  comes  in  analysis  of  each  window's  Fourier 
transform. 

If  a  window  is  chosen  which  has  a  very  narrow  main  lobe  and 
virtually  no  side  lobes,  then  the  data  can  be  passed  through  that 
window  and  the  resultant  frequency  transform  will  not  be  smeared.  One 
often  used  window  is  a  Hanning  window.  The  effect  of  windowing  is 
shown  in  Figure  9. 

Technically  the  data  is  now  distorted  and  the  effect  would 
corrupt  the  transfer  function.  However,  since  both  the  input  and 


output  data  signal  must  be  "windowed",  this  distortion  of  windowing  is 
removed  in  the  calculation  of  the  transfer  function  as  long  as  the 
same  window  is  applied  to  both  input  and  output  data. 

Figure  10  shows  various  windows  and  their  respective  Fourier 
transform.  The  question  of  which  window  is  best  for  this  application 
is  difficult  to  answer  since  several  factors  are  involved.  The 
"specifications"  of  the  Hanning  window  are  that  the  first  lobe  is  23 
dB  down  from  the  main  lobe  and  each  adjacent  lobe  falls  off  at  -12  dB 
per  octave  (40  dB  per  decade). 

Because  of  the  presence  of  these  sidelobes  there  is  still  some 
leakage  as  occurred  with  the  rectangular  window.  The  specifications  of 
the  original  rectangular  window  are  that  the  highest  sidelobe  is  13  dB 
down  from  the  main  lobe  and  each  sidelobe  falls  off  (decays)  at  6  dB 
per  octave.  Thus  the  Hanning  window  is  better  than  the  rectangular 
window  and  one  will  acheive  much  more  accurate  frequency  information 
if  the  input  and  output  signals  are  passed  through  a  window  such  as  a 
Hanning  window  prior  to  the  calculation  of  Fourier  transform. 

One  desires  a  window  with  a  very  large  dB  drop  between  the  main  lobe 
and  the  first  side  lobe  where  the  most  leakage  will  take  place.  In  addition 
each  successive  lobe  should  have  a  healthy  rate  of  dB  drop  from  the 
previous  lobe.  Otherwise  there  would  be  a  minute  leakage  from  signals 
far  removed  from  that  in  questioned  and  this  minute  leakage  would 
accumulate  and  distort  the  results. 

The  main  lobe  is  also  an  important  factor  in  selecting  the 
optimum  window.  One  desires  a  window  with  a  very  narrow  main  lobe. 

It  should  be  no  wider  than  the  frequency  resolution  required  for  the  aircraft 
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analysis.  Table  1  contains  a  listing  of  various  windows  and 
specifications.  The  effects  of  windowing  are  always  present,  and  one 
can  almost  always  acheive  better  results  by  using  a  window  other  than 
the  rectangular  window,  (i.e.  no  windowing) 

Window  Parameters 

The  characteristics  of  a  particular  weighting  scheme,  or  window, 
depends  in  general  on  the  number  of  data  points  used.  As  a  comparison 
of  windows,  (Ref  7),  Table  1  contains  a  list  of  some  pertinent 
parameters  based  upon  a  sample  size  of  50  data  points  (N-50).  The  3.0 
and  6.0  dB  bandwidth  is  indicated  by  frequency  bin  widths  and  indicate 
the  points  at  which  the  gain  of  the  shaped  filter  (window)  is  down  3 
or  6  dB.  The  overlap  correlation  is  based  on  a  correlation 
computation  given  in  Ref  7  and  is  a  relative  measure  of  the  degree  of 
correlation  between  random  components  in  successive  windows  as  a 
function  of  the  amount  of  overlap. 

Another  very  good  example  of  the  effect  of  windowing  data  is  to 
examine  a  signal  which  consists  of  only  two  frequencies.  One 
frequency  is  dominant,  the  other  is  40  dB  lower  in  signal  strength  but 
is  located  very  close  to  the  dominant  signal.  The  effects  of  various 
windows  are  shown  in  Figure  11. 

Thus,  the  type  of  window  of  filter  shape  chosen  for  the  frequency 
transform  of  the  time  data  has  a  significant  impact  on  the  detection 
of  a  weak  signal  located  close  to  a  strong  signal.  Good  detection  of 
these  weaker  signals  comes  from  using  windows  or  filters  which  have 
low  sidelobe  levels.  However  a  characteristic  of  these  filters  is 
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that  as  the  sidelobe  level  goes  down  (less  leakage)  the  mainlobe  gets 
wider  (poor  frequency  resolution). 

In  the  work  presented  in  here,  the  Blackman  window  was  used  to 
filter  the  input  and  output  signals.  The  Blackman  filter  is  easy  to 
implement  in  the  time  domain  (one  line  in  a  Fortran  program)  and  has 
very  good  specifications  in  terms  of  sidelobe,  sidelobe  fall-off, 
scalloping  loss,  etc. 

The  Hanning  window  is  typically  used  in  linear  system  testing; 
namely  vibrational  analysis.  The  reason  the  Hanning  window  is  popular 
is  that  its  specifications  are  generally  adequate  for  most 
applications  and,  it  is  very  easy  to  implement  in  either  the  time 
domain  or  the  frequency  domain.  Users  of  the  Hanning  window  find  it 
convenient  to  transform  input  and  output  signals  into  the  frequency 
domain  and  then  remove  the  leakage  effects  by  smoothing  the  frequency 
data  with  the  Hanning  filter.  Because  of  its  simplicity  and 
versatility  the  Hanning  filter  has  been  "hard-wired"  into  frequency 
analyzers-devices  specifically  design  to  perform  very  high  speed 
spectrum  analysis  and  identification.  It  should  be  noted  however, 
that  the  selection  of  the  "best"  window  or  filter  for  most  purposes 
assume  a  single  pass  of  the  entire  input  and  output  signals  to 
calculate  the  various  frequency  functions  Gx,  Gy,  Gxy,  Hxy,  etc.  In 
reality,  when  the  data  is  sectioned,  processed,  and  averaged  to  obtain 
results  then  the  apparent  advantage  of  one  window  over  another  may  be 
lost  in  the  averaging  (Ref  9). 


III.  FAST  FOURIER  TRANSFORMATION 


Introduction 

The  mathematical  foundation  used  in  the  frequency  analysis  of 
flight  test  data  is  the  Fourier  transform.  As  shall  be  shown,  the 
Fourier  transform  is  an  integral  part  in  the  calculation  of  input  and 
output  density  functions,  (Gx,  Gy),  cross  correlation  function  (Gxy), 
and  the  actual  frequency  response  function  (Hxy).  The  Fourier 
transform  can  be  used  for  characterizing  linear  systems  and  for 
identifying  individual  frequency  components  of  a  continuous  waveform, 
(see  Figure  3). 

A  very  good  interpretation  of  the  Fourier  Transform  is  given  by 
Brigham  (Ref  5).  The  Fourier  transform,  in  essence  separates  the 
waveform  into  a  number  of  sinusoids  at  their  respective  frequencies. 
The  sinusoids  can  be  summed  to  return  the  original  waveform.  The 
Fourier  transform  is  represented  by  the  sinusoids.  The  time  domain 
signal  gives  amplitude  and  time  information  about  the  signal  whereas 
the  frequency  domain  give  amplitude  and  frequency  information.  The 
Fourier  transform  is  given  by: 

X(«o)  -  J’x(t)e"j21r‘°tdt  (1) 

where  x(t)  ■  waveform  to  be  decomposed, 
and  X(to)  =  Fourier  transform  of  x(t) 

In  the  case  of  flight  test  data,  one  does  not  have  a  mathematical 
expression  representing  the  input  signal.  In  other  words,  no 
expression  exists  which  can  be  substituted  into  (1)  to  calculate  the 
Fourier  Transform  of  the  input  (or  output)  signals.  Instead  the  input 


and  output  consist  of  a  sequence  of  samples  of  measurment  transducers, 
taken  at  a  reqular  interval,  e.g.  every  .05  seconds.  Thus  instead  of 
a  continuous  signal  there  is  a  series  of  sampled  data  points.  If  a 
total  of  N  of  these  sample  data  points  are  obtained,  each  taken  at  T 
seconds  apart  (i.e.  T  =  sample  interval),  then  the  Discrete  Fourier 
transform  is  given  as: 

N-l 

X(n)  =  1/N  21  x(k)e"j2'Trnk/N  n=0,l  ,2,3 ,4 ,5 . N-l  (2) 

k=0 

The  resulting  Discrete  Fourier  transform,  X(n),  of  the  signal  x(k) 
will  contain  both  a  real  and  an  imaginary  part. 

An  N-point  series  in  the  time  domain  will  transform  into  an  N/2 
point  series  in  the  frequency  domain  as  the  last  N/2  points  will  be 
nothing  more  than  the  mirror  image  (the  minus  frequencies  values)  of 
the  first  N/2  points.  The  Digital  Fourier  transform  (DFT)  as  shown  in 
(2)  can  be  rewritten  in  the  form: 

N-l. 

X(n)  -  1/N  2Z  x(k)Wkn  where  W  =*  e  2  "  j/N  (3) 
k-0 

The  direct  calculation  of  (3)  require  the  solution  of  a  matrix  of 
terms, 

[X(n)J  -  [Wkn]  (x(k)] 

or  if  the  matrix  notation  is  expanded,  using  a  sample  size  of  four; 


X(0) 

111  1 

X(0) 

X(1) 

1  w1  w2  w3 

x(l) 

X(2) 

1  w2  w°  w2 

x(2) 

X(3) 

1  w3  w2  w1 

x(3) 

Fourier  Transform  Points  ■  matrix  made  of  x  time  series 

(i.e.  frequency  data)  e-j2lT/N  terms  data 
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The  matrix  of  e  J  *'  terms  makes  use  of  the  relationship  Wnfc  * 

Wn^m°d(^).  nkmod(N)  is  the  remainder  after  division  of  nk  by  N.  Thus 

if  N  =  4,  n=2,  and  k=3  then  =  W^.  (e- =e- ^  )  To  perform  the 

matrix  multiplication  shown  would  require  16  complex  products  and  16 

complex  additions.  Or  in  general  for  N  data  points,  N  complex 

multiplications  and  additions  would  be  required.  Cooley,  Turkey,  and 

others  (Ref  5)  have  shown  that  it  is  possible  to  make  efficient  use  of 

the  symmetry  present  in  the  above  matrix  to  reduce  the  number  of 

calculations  dramatically.  These  methods  of  computing  the  DFT  are 

called  Fast  Fourier  Transforms  (FFT).  Whereas  an  N  point  sequence 
2 

would  take  N  calculations  using  conventional  DFT's,  the  same  sequence 
could  be  calculated  using  FFT  techniques  with  only  (N/2)log2N 
computations.  Or; 

N  =*  1024  DFT  requires  1,048,576  complex  multiplications 

N*  1024  FFT  requires  5120  complex  multiplications 

Figure  12  reflects  these  results.  Thus  the  FFT  greatly  reduces  the 
amount  of  processing  time  required  to  transform  large  amounts  of  data. 

For  N»2®  ,  (  Y»  integer),  the  FFT  algorithms  are  simply  a 
procedure  for  factoring  an  N  x  N  matrix  into  Y  matrices  each  of  which 
are  also  N  x  N.  The  Y  matrices,  will  be  highly  symmetric  and 
contain  an  extremely  large  number  of  zeros.  FFT  algorithms  make  use 
of  this  symmetry,  the  large  number  of  zeros  in  the  matrices,  and  the 
fact  that  certain  patterns  are  even  present  for  the  data  in  "bit 
reversed”  order  to  reduce  the  number  of  calculations  involved.  A 
number  of  FFT  algorithms  are  available  for  performing  a  rapid  DFT.  A 
modification  of  an  algorithm  suggested  by  Brigham  (Ref  5)  is  employed 
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in  the  FFT  Fortran  subroutine  shown  in  Figure  13. 

Processing  Real  Data 

To  employ  the  subroutine  of  Figure  13,  it  should  be  recalled 
that  the  Fourier  transform  requires  a  complex  input  signal  and 
therefore  both  a  real  and  imaginary  portion  of  the  waveform  must  be 
provided.  Real  inputs  signal  should  have  their  corresponding  imaginary 
parts  set  equal  to  zero.  However  this  approach  is  wasteful  of  both 
computer  time  and  space.  In  the  actual  calculation  of  the  transforms 
of  input  and  output  time  series,  one  can  compute  the  FFT  of  both 
signals  simultaneously  by  inserting  one  signal  into  the  imaginary 
component  of  the  other,  transforming,  and  then  unscrambling  the 
results.  For  example,  if  x(k)  is  the  input  signal  and  y(k)  is  the 
output  signal,  where  the  input  and  output  signals  are  real  signals, 
then  one  normally  would  set  the  imaginary  part  equal  to  zero  and  then 
use  the  FFT  subroutine.  But  a  new  function,  z(k),  can  be  created  such 
that 

z(k)  =  x(k)  +  j  y(k) 

Through  the  use  of  the  linearity  property  of  the  Fourier 
transform  and  by  decomposing  frequency  signals  into  their  real  and 
imaginary  parts,  it  can  be  shown  (Ref  5)  that  the  input  and  output 
signal  can  be  recovered  in  the  transformed  state  by  the  following; 

X(n)  ■  (R(n) / 2  +  R(N-n)/2)  +  j  (IM(n)/2  -  IM(N-n)/2) 

Y(n)  -  (IM(n)/2  +  IM(N-n)/2)  -  j  (R(n)/2  -  R(N-n)/2)) 


where  Z(n)  ■  R(n)  +jIM(n) 

Thus  since  both  the  input  and  output  signals  are  real,  they  can  be 
transformed  at  the  same  time  with  virtually  a  50%  savings  in 


C  XXXXXXXXXXXXXXXXXXXXXXZXXXXXXXXXXXXXXXXXXXXXXXX 
C  SUBROUTINE  FFT  CALCULATES  THE  FOURIER  TRANSFORM 
C  OF  THE  INPUT  AND  RESPONSE  SIGNALS 
C 

SUBROUTINE  FFTCXREAL . XIMAG. N, NU. J ) 

DIMENSION  XREALCN).  XIMAGCN) 

C  ARRAY  XREAL  CONTAINS  REAL  PART  OF  TRANSFORM 
C  ARRAY  XIMAG  CONTAINS  IMAGINARY  PART  OF  TRANSFORM 
C  INPUT  DATA  IS  DISTROYED 

C  N  IS  NUMBER  OF  DATA  POINTS  (INTEGER  POWER  OF  2) 

C  NU  IS  POWER  OF  2  SUCH  THAT  N-2XXNU 

C  J  •  +1  GIUES  FORWARD  TRANSFORM  (-1  GIUES  INUERSE) 
N2-N/2 
NU1-NU-1 
K-0 

DO  100  L-l.NU 

102  DO  101  I-1.N2 

P* I B I TR( K/2XXNU1 • NU ' 

ARG-6 . 293 1 85  X  P / FLOAT C  N i 

C-CCSCARG5 

S-SIN(ARG) 

Kl-K+1 

K1N2-K1+N2 

TREAL  •  XREAL  C  K 1 N2 )  *C +X I  MAG  (  K 1 N2 )  *S 
7 ; nAG-x I MAG( K1N2 ) *C-XREAL ( K1N2) *S 
XREAL ( K1N2 ) -XREAL ( K1 ) -TREAL 
XIMAG(K1N2)-XIMAG(K1 5-TIMAG 
XREAL ( K1 ) -XREAL  C  K1 5 +TREAL 
X I MAGC K1 ) -XIMAGC  KD+TIMAG 
K-K+l 

101  CONTINUE 
K-K+N2 

IFCK  .LT.  N)  GO  TO  102 
K-0 

NU1-NU1-1 
N2-N2/2 
100  CONTINUE 

DO  103  K-l.N 
I-IBITRCK-l.NUJ+l 
IF  (I  .LE.  K)  GO  TO  103 
TREAL -XREAL ( K ) 

TIMAG-XIMAGCK) 

XREAL  ( K )  -XREAL  ( I ) 

XIMAGC K) -XIMAGC I ) 

XREAL ( I ) -TREAL 
XIMPG( I ) -TIMPG 

103  CONTINUE 

IF  (J  .EQ.  -1)  J-N 
DO  104  I-l.N 
XREAL ( I ) -XREAL ( I )/J 
XIMAGC I ) -XIMAGC I)/J 

104  CONTINUE 
RETURN 
END 

C 

C  FUNCTION  IBITR  PERFORMS  BIT  REUERSING  FUNCTION 
C 

FUNCTION  IBITRCJ.NU) 

Jl-J 
IBITR-0 
DO  200  I-l.NU 
J2-J1/2 

IBITR- IBITRX2+C J1-2XJ2 ) 

J1-J2 

200  CONTINUE 
RETURN 

END 


II  r'T  0  r  r-  «■  r  n 
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IV.  FREQUENCY  ANALYSIS  PARAMETERS 


Introduction 

To  analyze  a  system  in  the  frequency  domain,  several  parameters 
are  required.  One  needs  to  generate  the  Power  Spectral  Density  (PSD) 
of  the  input  and  output  signals.  The  PSD  of  the  input,  Gx,  is 
important  as  it  allows  one  to  examine  the  frequency  content  (spectrum) 
of  input  signal  to  tell  at  which  frequencies  the  input  signal  contains 
the  bulk  of  its  power.  For  proper  frequency  response  analysis,  it  is 
important  that  the  system  in  question  be  excited  across  a  wide  enough 
bandwidth.  Likewise  a  plot  of  the  spectrum  of  the  output  signal,  Gy, 
or  PSD  of  the  output  would  show  the  magnitude  of  the  frequencies  found 
in  the  output. 

Quantities  to  be  calculated  in  the  frequency  response  analysis 

are: 

1.  The  Power  Spectral  Density  (PSD)  of  the  input  signal,  x(t).  The 
PSD  of  the  input,  annotated  Gx,  is  sometimes  called  auto  spectral 
density  function  or  normalized  power  spectrum. 

2.  The  Power  Spectral  Density  of  the  output  signal,  y(t),  is 
annotated  as  Gy. 

3.  The  Cross  Spectral  Density  function  (CSD)  relates  the  input  and 
output  time  signals  together:  i.e.  the  output  function  y(t)  and  the 
forcing  function  x(t).  The  CSD  of  these  two  signals  is  denoted  as 
Gxy,  and  the  value  of  Gxy  is  used  in  the  calculation  of  the  Frequency 
Response  Function  and  the  coherence  function. 

4.  The  Frequency  Response  Function,  or  often  referred  to  as  the 
transfer  function,  will  be  denoted  by  Hxy. 


5.  The  coherence  function,  9  ,  will  be  calculated  and  is  helpful  in 
analyzing  systems  in  the  presence  of  noise. 

Armed  with  the  five  pieces  of  information  listed  above,  Gx,  Gy, 
Gxy,  Hxy,  and  ^  ,  one  can  now  provide  a  detailed  analysis  of  a  system 
(i.e.  aircraft,  flight  controls,  etc).  Below  is  a  description  of  each 
of  the  terms  and  how  they  are  generated  from  flight  test  data  for 
this  project. 


POWER  SPECTRAL  DENSITY  FUNCTIONS  (Gx  and  Gy) 

The  input  power  spectral  density,  Gx,  is  the  magnitude  squared 
of  the  complex  coefficients  of  the  Fourier  transform  of  x(t). 
Similarly,  Gy,  is  the  magnitude  squared  of  the  complex  coefficients  of 
the  Fourier  Transform  of  y(t).  For  example,  if  the  Fourier  Transform 
of  x(t),  evaluated  at  a  frequency  -2  n/N  is  written  as: 

xCWj)  =  +  j  bL 

then  Gx(u»^)  *  X(wi)X*(*Ji) 
where  X*(wi)  *  aA  -  jbj^ 

Therefore  GxCtj^)  -(a^)2  +(b^)2  (4) 

Similiarly,  if 

Y(*i)  -  <4  +  jdt 

then  Gy(*»>^)  ■  (Cj)2  +  (d^.)2  (5) 

The  power  spectral  density  is  an  easy  way  to  examine  a  waveform  to 
determine  what  frequencies  are  present  and  in  what  magnitudes.  Since 
the  analysis  of  aircraft  flight  control  systems  are  generally  most 
concerned  with  the  frequency  range  of  0.1  to  10  rad /sec,  one  needs 
to  examine  the  PSD  of  the  forcing  function  x(t).  Satisfactory  results 
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may  seemingly  be  obtained  from  an  unsatisfactory  flight  control  system 
but  close  examination  may  show  the  flight  control  system  was  never 
excited  in  the  frequency  range  where  the  problem  exists.  If  one 
expects  a  problem  near  say,  9  rad/sec  due  to  insufficient  damping,  and 
yet  the  input  shows  very  little  frequency  content  near  9  rad/sec 
occurred,  then  the  pilot  and  the  engineer  might  otherwise  deduce  that 
control  is  adequate  in  that  frequency  region. 

Cross  Spectral  Density  Function 

The  Cross  Spectral  Density  function,  Gxy,  requires  the  Fourier 
transform  of  the  input  x,  and  output,  y.  As  before,  if 
X  (cOi )  =  ai  +  jb^  and 

Y Ca»±)  *  cA  +  jd^ 

then  Gxy  is  defined  as  the  cross  conjugate  product  or 

GxyCtop  -  (at  -  jbi)(ci  +  jdt) 

expanding:  Gxy(tt»i)  »  (a^  +  bidi)  +  j(aidi  ~  bici) 

In  the  software  designed  for  this  project,  the  real  part  of  Gxy  is 
designated  P  and  the  imaginary  part  is  designated  Q.  Therefore; 

Pi  "  aici  +  bidi 
Ql  -  aid1  -  biCi 
or  GxyCw^  -  ?L  +  jQt 


Since  Gxy  is  a  complex  number,  it  contains  both  magnitude  and  phase 


information,  given  by: 


magnitude  of  Gxy^)  =  | Gxy ) | 
phase  of  GxyCt^)  =  ^^^cy(j^)  = 


(Pi2  +  Qi2)1/2 

tan^'tQj/Pi) 


Thus  one  can  calculate  the  magnitude  and  phase  of  Gxy  at  any 
frequency,  U)  The  importance  of  this  magnitude  and  phase 
information  is  shown  below. 


Frequency  Response  Function 

The  Frequency  Response  Function,  HxyCu^),  is  depicted  by  Figure 
14.  The  frequency  response  function  relates  the  input  to  the  output 
in  the  frequency  domain,  and  is  defined  as  HxyC*^)  =  Ydu^J/XO*^). 

The  frequency  response  function  is  simply  the  Fourier  transform  of  the 
impulse  response  function.  Note  that  by  using  the  previous 
definitions  of  Y  and  X: 

Hxy(u^)  =*  (ct  +  jd± )  /  (ai  +  jbi ) 

Since  Hxy  is  a  complex  quantity,  one  can  calculate  its  magnitude 
and  phase.  Multiplying  the  numerator  and  denominator  by  the  complex 
conjugate  of  the  denominator,  one  obtains: 


Hxy  -  (cj  +jdi)/(ai  +  jbi)  *  (at  -  jb1)/(ai  -  jbt) 


Hxy  *  [(a^  +  bid1)  +  jCa^  -  bici)]  /  (ai2  +  bj2) 


Therefore  the  phase  of  Hxy  is: 


/^xy  - 


tan-1  (a1d1  -  b1c1)/(aici  +  bid1) 


Note  that  this  is  exactly  the  same  expression  obtained  for  the  phase 
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of  Gxy.  The  magnitude  of  Hxy  can  be  found  as  follows: 


Magnitude  of  Hxy  -  (ci2  +  di2)1/2  /  (at2  + 

=■  (Gy)1/2  /  (Gx)1/2  (7) 

This  can  be  related  to  Gxy  by  expanding  Gxy  in  terms  of 
a  ^ ,  b ^ ,  c ^ ,  and  d ^ • 

Magnitude  of  Gxy  *  (P^2  +  Q^2)^2 

-  I(aiCl  +  b^)2  +(a1d1  -  biCi)2]1/2 

*  [(alCl)2  +  Za^c^  +  (b^)2  +(a1di)2 
-  2aibicidi  +  (biCl)2]1/2 

-  [(aiC;L)2  +  (b^)2  +  (a^)2  +  (blC;l)2]1/2 

-  [ai2(ci  +  dt)2  +  b.2(Ci  +  d^2]172 

-  ((3i2  +  bi2)(ci2  +  dt2)]1/2 

By  equations  (4)  and  (5),  the  product  terms  are  Gx  and  Gyt  and  therefore 
Magnitude  Gxyfc^)  -  [Gx(wi)Gy(fc>i)] 1//2 
dividing  both  sides  by  Gx: 

Magnitude  of  Gxy(Ui)/Gx(«ai)  ■  GyCta*^)1^2  /  CxC*^)1^2  (8) 

But  since  this  is  exactly  the  expression  for  the  Frequency  Response 
function,  Hxy,  as  given  in  (7),  Hxyfcp  -|Gxy(«i)|  / 1 GxC«»£>  | 

Thus  there  is  no  need  to  directly  calculate  the  frequency  response 


function  by  use  of  (7).  Instead  the  magnitude  of  the  frequency 
response  function  can  be  obtained  from  the  ratio  of  the  magnitude  of 
Gxy  and  Gx.  The  phase  angle  of  Hxy  is  simply  the  phase  angle  of  Gxy. 
The  processes  for  the  calculation  of  PSD,  CSD,  and  Frequency  Response 
functions  are  shown  in  Figures  15,  16,  and  17. 


NOISE  REDUCTION  TECHNIQUES 

The  major  source  of  error  in  the  frequency  response  plots  comes 
from  high  frequency  signals  corrupting  the  data  because  of  the 
aliasing  phenomenon  and  from  the  many  sources  of  noise.  Noise  is 
present  in  the  input  and  output  signal  measurement.  Noise  is  also 
introduced  in  the  digitizing  process  and  is  a  function  of  the 
resolution  of  the  measurement  sensors. 

For  this  program  the  coherence  function  was  used  to  provide 
information  to  the  user  at  which  frequencies  the  response  was  least 
corrupted  by  noise.  In  addition,  a  noise  reduction  algorithm  was  added 
to  the  basic  frequency  response  program  to  reduce  the  effects  of 
noise.  This  technique  and  the  coherence  function  are  described  below. 


With  the  transfer  function  measured  using  the  method  previously 
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discussed,  the  coherence  function  denoted  (  5  )  can  be  computed.  The 
coherence  function,  is  defined  as; 


U  _  (response  power  caused  by  applied  input) 
(measured  response  power) 


Ideally  in  the  absence  of  noise  the  coherence  function  would  be  unity. 
When  the  measure  response  power  is  greater  than  the  measured  input 
power,  e.g.  because  of  some  extraneous  noise  source  is  contributing  to 
the  output  power,  then  the  coherence  value  will  be  less  than  one. 


|"*3TA  !  N 

F  0  L'  K  1  E  K  COEFFICIE  X  T  S 
.i  AND  b  AT  FREQUENCY 


i  FROM  x  (  t  ) 


SQUARE  REAL  AND 
IMAGINARY  COEFFICIENTS 

COMPUTE  SUMMATION  OF 
SQUARED  TERMS 


I 

G  x  (to  .  )  =  a  .  +  b 
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This  will  be  the  case  for  those  frequencies  where  the  noise  source 
adds  power  to  the  response  signal. 

Therefore  the  coherence  function  can  be  used  to  indicate  the 
degree  of  noise  contamination  in  the  transfer  function  estimation. 
Because  flight  test  data  will  inherently  be  corrupted  by  noisy 
measurements,  the  evaluation  of  IS*’  for  the  frequency  range  0.1  to  10 
rad/sec  is  extremely  valuable.  Figure  18  shows  an  example  of  the 
noise  problem.  Since  one  is  most  interested  in  identifying  modal 
parameters  from  measured  transfer  functions,  the  variance  of  the 
estimates  of  those  parameters  is  reduced  in  proportion  to  the  amount 
of  noise  reduction  in  the  measurements. 

From  Figure  18,  one  can  compute  H(w)  from  the  measured  and 
transformed  signals  X(w)  and  M(««»).  H(«o)  is  defined  as  the  ratio 
Y(u)/X(cJ).  Thus; 

M(«j)  =  H(oa)X(o>)  +  N(u>)  or  simply  M=HX  +  N 

The  input  power  spectrum  is  Gxx  and  is  defined  as  Gxx  =  XX* 
where  denotes  the  complex  conjugate  of  the  transform.  The  cross 
power  spectrum  is  Gmx  *  MX*  *  (HX  +  N)X  =  HGxx  +  NX* 

Consider  averaging  the  quantity  Gmx.  The  average  value  of  Gmx  from  n 

n 

different  measurements  is  Gmx  *  l/n^Gmx(i)  where  Gmx(i)  is  the  ith 

1  __ 

measurement  taken.  Therefore  Gmx  -  HGxx  +  NX  =  HGxx  +  Gnx. 

Likewise  Gxx  and  Gnx  are  defined  in  a  similar  manner.  Therefore  the 
actual  frequency  response  of  the  system  is  given  by; 

Hmx  *  (Gmx/Gxx)  -  (Gnx/Gxx) 

notice  that  as  the  number  of  averages  grows  larger  the  noise  term  is 
reduced  and  the  ratio  Gmx/Gxx  more  accurately  estimates  the  true 


Fourier  Transformation  of  Input  Signal 
Frequency  Response  Function 
Fourier  Transform  of  Output  Signal 
Noise 

Fourier  Transform  of  Measured  Signal 


Figure  18  Measurement  of  Signal  Plus  Noise 


transfer  function.  A  similar  analysis  is  presented  in  Ref  6  and  12 
for  the  frequency  response  of  structures. 

Figure  19  shows  a  frequency  response  of  roll  rate  to  stick  force 
for  a  YA-7D  (Ref  10).  Figure  19  contains  no  averaging  techniques  and 
notice  how  at  the  higher  frequencies  much  uncertainty  exists  as  to  the 
exact  nature  of  the  magnitude  and  phase  information.  Notice 
especially  how  the  phase  information  is  adversely  affected  by  noise. 
Figure  20  is  a  frequency  response  plot  of  stick  force  to  roll  rate 
using  the  same  data  as  was  used  to  generate  Figure  19.  The  program 
used  to  generate  Figure  20  used  the  same  basic  techniques  as  the 
program  used  to  generate  Figure  19,  however  instead  of  examining  a 
single  512  data  point  sequence,  the  program  divided  the  512  data 
points  into  seven  128  data  point  segments  as  shown  in  Figure  21. 

While  some  frequency  resolution  was  lost  using  this  technique, 
the  resolution  was  still  high  enough  for  aircraft  parameter 
determination.  More  importantly,  the  averaging  of  the  seven  segments 
of  data  significantly  reduced  the  scatter  of  the  magnitude  and  phase 
curves  especially  in  the  lower  range  of  frequency.  This  technique  has 
now  provided  usable  magnitude  and  phase  information  in  the  frequency 
range  of  interest,  namely  0.1  to  10  radians  per  second.  Larger  sample 
sizes  with  even  more  segments  would  reduce  noise  even  further.  Notice 
that  even  with  average  techniques  the  curves  are  not  usable  beyond 
about  15  radians  per  second.  This  is  due  to  the  aliasing  problems 
described  earlier.  It  is  also  a  function  of  the  sampling  rate. 
Engineering  experience  shows  that  reliable  frequency  information  is 
not  achieved  at  the  Nyquist  limit  of  2  samples  per  cycle.  Useful 
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Division  of  a  512  Point  Data  String  into 
Seven  128  Point  Data  Segments 


information  may  not  be  obtained  in  some  cases  until  at  least  five 
samples  per  cycle  occur.  In  the  case  of  a  0.05  second  sample  rate  (20 
samples  per  second),  the  maximum  usable  frequency  of  five  samples  per 
cycle  occurs  at  just  below  21  radians  per  second.  Figure  22  is  a  plot 
of  the  Coherence  Function  C8V)  for  the  data  used  to  generate  Figure 
20.  The  Coherence  function  has  been  scaled  by  a  factor  of  20 
therefore  frequencies  having  a  gamma  squared  value  of  20  indicate  very 
good  correlation  between  input  and  output  (no  noise).  Figure  22  shows 
that  above  10  radians/second  the  gamma  squared  or  coherence  function 
begins  to  fall  significantly  indicating  poor  correlation  between  input 
and  output  (noise).  The  gamma  squared  function  is  determined  by: 

^  =  Mag  Gxyz  /  (Mag  Gx  *  Mag  Gy) 

Plotting  Programs 

To  make  a  quick  look  examination  of  the  input  or  response  time 
histories,  a  program  called  PLOT  was  written  which  displays  to  the 
Test  Pilot  School  user  on  the  terminal  the  data  string.  The  plot 
program  is  excellent  for  determining  the  quality  of  the  time  history 
(excessive  noise,  low  frequency  content,  wild  points,  etc) 

The  frequency  response  information  from  the  program  FRAN 
(FRequency  ANalyzer)  is  placed  into  files  called  B0DE.DAT  (magnitude 
and  phase  information),  GX.DAT  (power  spectral  density  of  input),  and 
GAMMA2.DAT  (  or  coherence  function).  This  information  is  then 
plotted  on  semi-log  graphs  using  a  modified  Test  Pilot  School  plotting 
routing  called  BDEPLT.  The  computer  processing  for  all  program  is 
done  on  the  Test  Pilot  School  PDP-11/34  computer.  Hard  copy  plots  are 
produced  on  a  Gould  electro-static  plotter. 
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SIMULATED  FLIGHT  TEST  DATA 


To  verify  the  FRAN  program  which  performs  the  frequency 
analysis,  it  was  desired  to  test  the  program  not  with  actual  flight 
test  data  but  rather  with  simulated  data.  Several  programs  were  used 
to  generate  the  simulated  flight  test  data.  Program  F410  is  for  an 
F-4  at  55,000  feet  and  mach  number  of  1.8  M.  Program  F45  is  for  an 
F-4  at  35,000,  0.6M,  and  F89  is  for  an  F89  at  20,000  feet  0.7M.  The 
programs  are  all  identical  except  for  the  respective  A  and  B  matrices. 
Each  program  solves  the  basic  equation: 

x(t)  -  Ax( t )  +  Bu( t ) 
where  x^  *  state  vector 

u  ■  control  vector 

a  solution  of  the  above  is  expressed  as: 

x(t)  -  eAtx(0)  +  /  eA(t-'r>Bu<'T)d'r 
0 

If  t  ■  kT  (discrete  time  intervals) 

kT 

x(kT)  -  eAkTx(0)  +  S  eA(kT  BuCTJd'i' 

0 

for  the  case  of  k*2, 

T  2T 

x(2T)  -  eA^Tx(0)  +  /  eA(2T"T)Bu(0)d^  +  f  eA(2T_r)Bu(T)dr 

0  0 

it  can  be  shown  that 
eAkTx(0)  -  eATx[(k-l)T] 

so  in  general  and  using  the  identity  matrix  I: 


x [ (k+l)T ]  =  eAi x(kT)  +  A_1(eAi  -  I)Bu(kT) 


at 

e‘  was  replaced  by  a  series  expansion  and  approximated  by  the  first 
four  terms: 

eAT  =  [I  +  AT  +  A2T2/2!  +  A3T3/3!  ] 
giving; 

to  OS 

X[(k+1)T]  =  [Y1  AkTk/k!  ]x(kT)  +[  £L  AkTk+1/(k+l)!  ]u(kT) 
k=*0  k=0 

which  was  the  basic  equation  solved  in  all  the  programs  which 
generated  the  simulated  flight  test  data.  The  program  source  code  is 
Fortran  and  the  program  F410  is  shown  in  Appendix  B.  Figure  23  shows 
the  flowchart  for  generating  the  simulated  flight  test  data. 

PROGRAM  FEATURES 

The  program  requests  the  amount  of  output  (time)  from  the  user. 
The  program  also  requests  the  sampling  rate  (i.e.  0.125  seconds  per 
sample).  Typically  output  time  and  sampling  rates  were  chosen  to 
provide  a  convenient  amount  of  data  for  the  frequency  response  program 
(FRAN)  (i.e.  512  or  1024  data  points). 

While  the  program  generates  output  at  the  sample  rate  requested 
by  the  user,  it's  step  size,  T,  between  output  can  be  made  smaller  at 
the  discretion  of  the  user.  The  value  of  T  was  typically  0.01  seconds 
The  program  was  designed  with  A  and  B  matrices  representing  the  pitch 
axis.  The  output  of  this  program  was  the  change  from  equilibrim  of 
pitch  rate  (q),  pitch  angle  (©),  velocity  (u),  and  angle  of  attack, 

(0t  ).  An  example  of  the  output  is  shown  in  Figure  24.  Data  files 
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(i.e.  VEL0CITY.DAT)  are  created  for  each  of  the  program's  outputs. 

PROGRAM  FRAN  VERIFICATION 

To  verify  the  frequency  response  program  (FRAN)  the  program  was 

run  with  the  simulated  flight  test  data  generated  from  program  F410. 

(F-4  at  55,000  feet  and  1.8M).  The  F-4  at  this  flight  condition  is 

known  to  have  a  short  period  natural  frequency  (co  )  of  4.84  rad/sec 

sp 

and  a  short  period  damping  ratio  (*fSp)  of  0.06.  The  input  for  the 
F410  program  is  elevator  deflection  (  i  e).  A  "random”  input  signal 
was  created  by  arbitrarily  hand  drawing  and  digitizing  an  input 
signal . 

Simulating  a  system  which  provides  output  eight  times  per  second, 
a  sampling  rate  of  0.125  second  was  entered  along  with  a  request  for  64 
seconds  of  output.  This  combination  generated  512  data  points  which 
is  an  integer  exponential  of  2  (a  requirement  of  the  FRAN  program). 

The  input  signal  (512  points  of  e)  along  with  a  simultaneous  response 
signal  (512  data  points  of  velocity  changes  (u),  were  used  as  the 
signals  to  be  frequency  analyzed  by  the  program  FRAN.  The  magnitude 
and  phase  information  created  by  the  program  FRAN  were  plotted.  This 
plot  is  shown  in  Figure  25.  Note  that  Figure  25  shows  a  lightly 
damped  short  period  to  exist  near  4.3  rad/sec  and  the  DB  increase  at 
this  frequency  of  10  DB  along  with  the  180  degree  phase  shift, 
indicates  an  equivalent  second  order  damping  ratio  of  about  0.06. 

This  identification  of  a  short  period  mode  by  frequency  and  damping 
ratio  is  in  very  close  agreement  with  the  actual  values  indicating  the 
frequency  response  program  FRAN  did  in  fact  work  properly.  These 
values  were  also  in  close  agreement  with  the  eigenvalues  determined 
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FREQUENCY  -  RRD/SEC 


from  Che  A  matrix  which  generated  the  original  simulated  flight  test 
data.  The  improper  increase  in  the  curve  at  the  higher  frequencies 
(above  10  rad/sec)  appears  to  be  a  function  of  poor  phase  information 
due  to  so  few  sample  per  cycle  and  a  function  by  which  the  response 
signal  is  generated  in  the  program  F410  (not  truly  a  simultaneous 
input  and  response  signals).  This  phenomenon  was  confirmed  by 
increasing  the  sampling  rate  to  20  samples  per  second.  In  this  case 
the  phase  increase  was  still  present  but  a  a  much  higher  frequency- 
outside  the  frequency  range  of  interest  of  0.1  to  10  radians/second. 
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V.  LOWER  ORDER  EQUIVALENT  SYSTEM  DETERMINATION 


Once  the  frequency  response  plot  were  determined  using  the 
frequency  analysis  methods  of  the  program  FRAN,  the  next  step  was  to 
develop  a  means  of  reducing  the  complex,  nth  order,  pilot-in-the-loop 
aircraft  responses  to  a  lower  order  equivalent  system  (LOES). 

Background 

Lower  order  equivalent  systems  allow  the  engineer  to  compare  the 
responses  of  a  mth-order  aircraft  to  a  n-th  order  aircraft  by  reducing 
both  systems  to  a  common  order.  For  example,  responses  of  the  pitch 
axis  would  likely  be  reduced  to  4th  order  (phugoid  and  short  period 
modes)  or  even  just  second  order  (short  period  only).  Responses  in 
the  roll  axis  would  likely  be  reduced  to  a  4th  order  (spiral,  dutch 
roll,  and  roll  mode)  system  and  may  be  reduced  as  low  as  a  first  order 
(roll  mode  only).  Lower  order  equivalent  system  also  are  becoming  the 
basis  by  which  aircraft  will  be  compared  to  the  requirements  of  the 
"new"  MIL-SPEC  (Ref  8).  The  program  which  was  developed  for  the 
United  States  Air  Force  Test  Pilot  School,  provides  a  Last  fit"  of 
the  actual  flight  test  frequency  response  data  to  virtually  any 
transfer  function  the  user  may  desire.  For  the  pitch  axis,  the  most 
common  LOES  used  was 

0  K.q  [  s  +  (  1/Tq)  ]  e-Tds 

e  s  +  2^gpU^pS  +  u»Sp 
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response  and  response  produced  by  the  LOES.  A  subroutine  calculates  a 
gradient  and  a  gradient  search  begins  to  minimize  the  error.  A 
scaling  rule  is  incorporated  to  improve  search  convergence,  and  a 
constraint  is  placed  in  the  search  to  insure  model  stability  during  the 
search. 

The  software  which  was  developed  for  this  lower  order  equivalent 
system  fit  utilized  subroutines  developed  at  tne  Lewis  Research  Center 
by  Robert  C.  Seidal  (Ref  13).  The  subroutines  were  designed  to  be 
used  with  a  program  which  was  designed  to  fit  various  transfer 
function  models  to  given  experimental  frequency  response  data.  Because 
of  the  sophistication  of  the  subroutines  and  gradient  calculation,  it 
was  determined  that  use  of  the  subroutines  might  provide  a  good  means 
of  determining  lower  order  equivalent  system  fits  for  flight  test 
data . 

Program  Formulation 

The  lower  order  transfer  function  provided  by  the  user  is  matched 
to  actual  frequency  response  over  a  range  of  frequencies  provided  by 
the  user  and  is  typically  in  the  range  of  0.1  to  10  rad/sec.  The 
actual  frequency  response  Is  expressed  as  a  complex  number  A(«»>)e^  ^ 
where  A  is  the  magnitude  and  ©  is  the  phase.  The  lower  order 
equivalent  model  transfer  function  G(ju>,b),  where  b  is  a  vector  of 
parameters  supplied  initially  by  the  user.  The  cost  function  J(b)  is 
defined  as 

CD 

•1(b)  ■  /  [  G(jw,b)  -  A ((•»)  e-i*  ^  do>  (fc»  ■  frequency) 

0 

To  make  the  above  equation  suitable  for  the  computer  the  following 
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where  <ogp  =  equivalent  short  period  natural  frequency 
^sp  *  equivalent  short  period  damping  ratio 
=  equivalent  time  delay 
Kq  -  gain 


For  the  lateral-directional  axis  the  most  common  LOES  system  was 

JL  ,  KP  (0)  [  $o»“qI  e"TdS 
Fs  (l/Ts)(l/Tr)I  J  dr.  «OdrJ 


where  Tr  =  equivalent  roll  mode  time  constant  P  =  roll  rate 

Wdr  *  equivalent  dutch  roll  natural  frequency  Fg  =  stick  force 

^dr  =  equivalent  dutch  roll  damping  ratio  K  =  roll  rate  gain 

Td  =•  equivalent  time  delay  Tg  =  spiral  time  constant 

[  $  ]  *  damping  ratio  and  natural  frequency  of  quadratic  term 

For  many  applications  in  the  lateral  axis,  such  as  the  DIGITAC  HAVE 
DELAY  project  (Ref  10),  the  LOES  was  simplified  to  a  first  order 
system  with  time  delay  which  provided  good  results.  The  LOES  used 
was ; 

P  K  e“Tds 

Fs  (s  +  1/Tr) 


where  K  *  gain 

Td  *  equivalent  time  delay 

1  /'fr  *  pole  location  corresponding  to  the  roll  mode  time 
constant 

The  method  of  determining  the  best  fit  of  a  equivalent  system 
transfer  function  to  frequency  response  data  uses  a  program  which 
determines  a  cost  function  which  is  the  error  between  actual  frequency 
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approximation  was  used 


J(b)  -  I/25:  [GCj^.b)  -  A(«1)e-)®(“>i>]2  («i+1 
i-1 

where  Nd  =  number  of  data  points. 

The  choice  of  the  frequencies  for  the  data  is  important.  If  the 
frequencies  are  improperly  chosen,  the  value  of  the  error  at  a 
particular  frequency  could  have  a  disproportionate  effect.  In  certain 
cases  it  is  desired  to  favor  a  close  fit  between  the  model  and  the 
plant  over  a  portion  of  the  frequency  range.  For  handling  qualities 
evaluations  where  this  program  is  most  often  used,  the  range  of 
frequencies  used  was  either  0.1  to  10  radians  per  second  or  1  to  10 
radians  per  sec.  Probably  the  best  spacing  method  is  to  use  evenly 
spaced  frequencies  from  a  logrithmic  scale. 


g2(s,b)  -  b2  gain  (11) 
g3(s,b)  -  (s/b3k  +  1  )-1  pole  (12) 
g^s.b)  -  [(s2/b3k)2  +  2sb^/bj  +  1]  quadratic  zero  (13) 
g6(s,b)  -  [s2/b7k  +  2sb6k/b7k  +  1]_1 


quadratic  pole  (14) 


g8( s , b)  =  e~sb8  time  delay  (15) 

The  factors  of  equation  (10),  (12),  (13),  and  (14)  are  repeatable  in 
equation  (9).  The  elements  of  the  b  vector  are  the  zeros  blk,  the 
gain  b2,  the  poles  b3^,  the  quadratic  zero  damping  ratios  bA^,  the 
quadratic  zero  natural  frequencies  bS^,  the  quadratic  pole  damping 
ratios  bb^,  the  quadratic  pole  natural  frequencies  b7^,  and  the  lead 
or  lag  time  b8.  The  number  of  gm(s,b)  factors  in  the  transfer 
function  is  the  number  of  parameters  in  the  b  vector  (np)  minus  the 
number  of  quadratic  factors  (nq). 

GRADIENT  CALCULATION 

The  gradient  of  J(b)  is  required  by  the  gradient  search  routine 
to  find  the  parameters  of  b  which  minimize  J(b).  If  parameters  are 
not  scaled,  many  searches  will  not  converge.  Gradient  search 
techniques  generally  require  parameter  scaling  to  obtain  efficient 
search  convergence.  For  this  program,  a  scaling  law  was  incorporated 
such  that  parameters  tend  to  have  equal  influence.  Each  parameter  was 
scaled  by  its  own  magnitude.  The  gradient  calculation  and  scaling 
routines  were  provided  in  a  subroutine  by  NASA.  Large  changes  in  the 
certain  values  of  the  b  vector  during  the  search  may  cause  a  jump 
across  a  relatively  narrow  error  boundary.  The  subroutine  prevents 
values  from  crossing  zero  by  constraining  values  as  they  approach  zero. 
Parameters  about  to  cross  zero  are  assigned  a  very  small  value  on  the 
same  3ide  of  zero.  This  prevents  changes  between  right  and  left  hand 
plane  zeros,  and  it  prevents  changes  between  lead  and  lag  time 
exponentials . 
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PROGRAM  MECHANIZATION 

A  complete  listing  of  this  program  is  contained  in  Appendix  B. 

The  main  program  asks  the  user  if  he  wishes  to  use  hand  entered  data 
or  data  from  a  file  called  FREQ.DAT.  Should  the  user  desire  to  hand 
enter  data,  the  file  FREQ.DAT  is  automatically  created  for  later  use. 
The  program  then  asks  the  user  to  input  the  LOES  transfer  function 
which  is  to  be  used  for  the  frequency  matching.  Up  to  15  parameters 
in  the  model  are  available.  The  user  inputs  starting  values  for  each 
of  the  parameters  in  the  model.  Once  the  LOES  transfer  function  and 
the  starting  values  have  been  entered,  the  user  is  offered  a  menu  of 
programming  instructions.  The  menu  is:  INSTRUCTIONS?  1=SEARCH 
2= NEW  FREQUENCY  DATA  3=  NEW  TRANSFER  FUNCTION  .  4=  PRINT  AND  SAVE 
RESULTS  5  *  STOP.  Instruction  1  causes  the  search  routine  to  begin 
and  the  best  fit  LOES  is  determined.  Instruction  2  allows  the  user  to 
re-enter  or  change  frequency  data.  Hand  entered  data  is  inputted  by  a 
frequency  list,  and  the  corresponding  amplitude  and  phase.  Amplitude 
is  given  in  dB's  and  phase  is  in  degrees.  Instruction  3  reinitializes 
the  program  to  the  point  of  entering  the  transfer  function. 

Instruction  4  prints  out  the  frequency  response  of  the  LOES  after  the 
best  fit  has  been  accomplished  and  also  stores  the  results  in  a  file 
called  "RESULTS.DAT".  Instruction  5  stops  the  program. 

The  program  output  contains  the  value  of  the  cost  function  (F), 
the  number  of  times  the  cost  function  was  calculated,  (KNT),  the  final 
parameter  step  size  (size),  and  the  program  also  contains  information 
as  to  whether  or  not  the  search  converged  or  encountered  an  error. 
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Chapter  VI  contains  examples  of  the  use  of  the  program  LOES  with 
the  Frequency  Response  Data  generated  from  flight  test. 


VI.  RESULTS  OF  FLIGHT  TEST 


Introduction 

This  chapter  presents  the  partial  results  of  a  flight  test  project 
which  incorporated  the  frequency  response  methods  and  Lower  Order 
Equivalent  System  fit  methods  presented  in  this  paper. 

A  project  called  HAVE  DELAY  (Ref  10),  was  conducted  by  selected 
members  of  USAF  Test  Pilot  School  Class  84-B.  The  objective  of  the 
Have  Delay  Project  was  to  investigate  the  combined  effects  of 
equivalent  time  delay  and  roll  mode  time  constant  on  the  lateral 
handling  qualities  of  fighter  aircraft.  The  test  aircraft  was  a  YA-7D 
DIGITAC  which  was  a  modified  A-7D.  The  major  modification  to  the 
DIGITAC  aircraft  from  the  basic  A-7D  is  that  the  analog  control 
augmentation  system  (CAS)  was  replaced  by  a  digital  control 
augmentation  system  (DCAS).  The  aircraft  was  also  modified  with  a 
Yaw,  Angle  of  Attack,  Pilot-Statics  System  (YAPS)  head  mounted  on  a 
flight  test  nose  boom,  and  sensitive  flight  test  instruments.  The 
testing  was  conducted  in  the  cruise  configuration  with  6  MAU-12B/A 
pylon  racks,  an  instrumentation  pod  mounted  on  the  right  inboard 
station  and  an  Inert  MK-82  mounted  on  the  left  inboard  station  to 
maintain  weight  and  drag  symmetry. 

The  test  was  conducted  at  15,000  feet  (pressure  altitude)  and 
involved  tracking  a  target  aircraft  which  was  maneuvering  at  2.5  g  and 
350  KIAS.  The  test  was  conducted  in  accordance  with  the  limitation 
specified  in  the  HAVE  DELAY  Test  Plan,  the  A-7D  Flight  Manual,  and 
AFFTC  Regulation  55-2  (Ref  11,  2,  and  1). 

Twenty  test  sorties  totaling  30.5  hours  were  flown  between  2 


April  1985  and  8  May  1985  at  the  Air  Force  Flight  Test  Center  (AFFTC), 


Edwards  AFB ,  California. 

The  test  objectives  were  designed  to  evaluate  fighter  type 
aircraft.  Specifically  the  objectives  were: 

1.  Correlate  lateral  handling  qualities  ratings  (Cooper-Harper)  to 
various  changes  in  equivalent  time  delay  in  the  lateral  axis. 

2.  Correlate  pilot  handling  qualities  ratings  to  changes  in 
equivalent  roll  mode  time  constant. 

3.  Correlate  pilot  ratings  with  equivalent  bandwidth. 

To  accomplish  the  objectives  of  the  Have  Delay  test  program  a 
method  of  determining  the  Equivalent  Parameters  was  needed  as  the 
Test  Pilot  School  had  no  means  of  performing  frequency  analysis  on 
flight  test  data.  The  method  used  was  to  incorporate  the  frequency 
response  software  (FRAN)  and  the  Lower  Order  Equivalent  System  fit 
software  (LOES)  previously  described  and  develop  for  this  thesis. 

Test  Article  Description 

The  A-7D  is  a  single  engine,  single  place,  transonic  light 
surface  attack  aircraft  manufactured  by  the  Vought  Aeronautics 
Company.  The  A-7D  is  powered  by  the  Allision  TF41-A-1  non-af terburing 
turbofan  engine.  A  complete  aircraft  description  is  contained  in  the 
Flight  Manual  (Ref  2)  and  the  Partial  Flight  Manual  (Ref  4). 

Test  Instrumentation  and  Data  Reduction 

Flight  test  data  was  recorded  by  an  on-board  magnetic  tape. 
Flight  test  data  was  also  displayed  real-time  on  ground  based  strip 
charts  via  a  Telemetry  system  (TM).  Stick  force  and  roll  rate  time 


histories  were  examined  for  noise  and  wild  points  and  if  found 
suitable  were  then  used  as  the  input  and  response  signals  for  the 
frequency  response  analyzer  program  (FRAN). 

Test  Methods  and  Conditions 

An  A-7  aircraft  was  used  to  provide  a  maneuvering  target  for  the 
DIGITAC  aircraft  to  track.  The  maneuvering  performed  by  the  target 
was  a  "canned"  maneuver  and  did  not  vary  from  test  point  to  test 
point.  Project  pilots  were  given  an  unknown  combination  of  lateral 
time  delay  and  roll  mode  time  constant.  Each  pilot  then  aggressively 
tracked  the  maneuvering  target  and  provided  a  handling  qualities 
rating  for  that  configuration. 

The  frequency  analyzer  program,  FRAN,  provided  a  means  of 
examining  the  frequency  content  of  the  pilot's  stick  force  (PSD  of  Fs) 
to  ensure  the  target  was  being  aggressively  tracked  and  to  ensure  the 
aircraft  dynamics  were  being  adequately  excited  throughout  the 
frequency  range  of  interest  (out  to  10  rad/sec).  Figure  26  shows  the 
power  spectral  density  of  the  input  signal,  stick  force  for  two  such 
configurations.  The  first  configuration  being  tested  in  Figure  26  is 
for  an  equivalent  roll  mode  time  constant  of  0.25  seconds  and  a 
equivalent  time  delay  of  150  ms*  The  other  had  a  0.50  second  time 
constant  and  a  time  delay  of  75  ms. 

Figure  27  shows  the  frequency  response  information  (BODE  plot) 
for  roll  rate  to  stick  force  for  configuration  HD-000-050.  In  all,  27 
different  test  configurations  were  evaluated  (i.e.  27  different 
combination  of  roll  mode  time  constant  and  time  delay).  To  analyze 
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Figure  26  Power  Spectral  Density  of  Pilot's  Input 
Configurations  HD-000-020  and  HD-000-040 


the  results,  27  different  frequency  response  plots  were  generated. 

From  the  data  created  to  make  these  plots,  the  lower  order  equivalent 
system  fits  were  made. 

Examination  of  Results 

In  producing  accurate  lower  order  equivalent  system  fits  to 
actual  aircraft  response,  two  main  requirements  must  be  met.  First 
the  frequency  response  results  must  accurately  represent  the  aircraft. 
Secondly,  the  LOES  fit  must  sufficiently  fit  the  flight  test  data. 

The  accuracy  of  the  frequency  response  methods  are  dependent  on  the 
length  of  the  time  history  used,  data  sampling  rates,  noise  reduction 
techniques,  window  filter  design,  and  the  amount  of  noise  present  in 
the  flight  test  signals,  (see  chapters  2  and  4). 

The  accuracy  of  the  LOES  fit  is  simply  a  matter  of  examining  the 
aircraft  magnitude  and  phase  information  and  comparing  to  the  magnitude 
and  phase  information  of  the  lower  order  system.  The  cost  function  is 
also  an  excellent  tool  for  determining  the  adequacy  of  the  fit.  For 
example  in  the  Have  Delay  test  program,  the  correlation  between  pilot 
ratings  and  equivalent  system  parameters  appeared  to  be  very  poorly 
defined.  However,  examination  of  the  LOES  fits  revealed  a  number  of 
fits  with  excessively  large  cost  functions.  These  poorly  fit 
parameters  were  removed  and  then  the  correlation  between  pilot  ratings 
and  equivalent  system  parameters  was  clear  and  distinct  (see  Ref  10). 
The  problem  appears  to  be  that  the  assumption  that  the  data  be  fit 
with  a  first  order  response  and  time  delay  for  the  roll  rate  to  stick 
force  transfer  function  was  not  valid  in  some  cases.  In  fact,  due  to 
the  manner  in  which  the  DIGITAC  software  was  designed,  some  particular 


test  configuration  were  significantly  higher  order  and  did  not  reduce 
well  to  a  first  order  system.  Figure  28  shows  and  example  of  both  the 
actual  flight  data  and  the  LOES  data  generated  when  the  flight  test 
data  was  reduce  to  a  first  order  system  with  time  delay.  Figure  29 
shows  the  raw  frequency  response  data  and  the  LOES  frequency  response 
data  as  a  comparison. 
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YA-70  OIGITRC  HAVE  DELAY  7  MAY  85  MISSION  18 
FREQUENCY  RESPONSE  LATERAL  TRACKING  TASK 
ROLL  RATE  TO  STICK  FORCE  20  SAMPLES  PER  SECOND 


FREQUENCY  -  RAD/SEC 

quency  Response  for  Conf igurai io 
ar  Order  Equivalent  System  Kit  ( 


■J.'  H-"  ' 


T 


r-,'  f'.  *  7  wj  v-j~i 


k*  ■ 

c’- 

I 


frequency 
-ad  sec' 


0.73630 
0.98200 
1 . 22700 
1 . 47600 
2.20900 
2 . 45400 
3. 19100 
4. 17200 
5. 154O0 
6. 13600 
8.09900 
10.06000 


amplitude 

(dB) 


21.45400 

20.99910 

20.70940 

20.50880 

19.98620 

19.72240 

18.97010 

17.64070 

16.61700 

13.95190 

11.44260 

7.91984 


phase 

(deg) 


-35.61970 

-43.55120 

-49.00160 

-52.91250 

-62.20370 

-66.52180 

-78.11060 

-96.30300 

-124.22000 

-155.02000 

-174.55200 

-216.28200 


Frequency  Response  Data  for  Configuration  HD-075-055 
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f  requency 

amplitude 

phase 

v rad .  sec  ) 

(dB) 

(deg) 

0.73630 

21.26440 

-23.08000 

0.98200 

21.06311 

-30.54047 

1.22700 

20.86127 

-37.77857 

1.47600 

20.60031 

-44.91162 

2.20900 

19.60699 

-64.50171 

2.45400 

19.35183 

-70.58494 

3. 19100 

18.31239 

-87.62411 

4.  1  ""*200 

16.94656 

-107.89097 

5.  15400 

15.66860 

-126.20632 

6.  13600 

14.50053 

-143. 14619 

0.09900 

12.4746 2 

-174.43457 

10.06000 

10.79648 

-203.66895 

Frequency  Response  Data  for  Lower  Q^der  Equiva)**-*  . 


K  e~Tds 

DDES  =  - 

(s/ ( i/Tr )  *  1) 


=  equivalent  roll  mode  time  constant 

Tj  =  equivalent  time  d»lav 

Cost  Function  =  27.55064 
lain  (K)  =  11.8972 
’^r  =  3.0598 
",  =  226  ms 


Figure  29  Frequency  Response  Data  and  Data  for  Corresponding 
Lower  Order  Equivalent  System  Fit  (LOLS) 
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VII.  RESULTS,  CONCLUSIONS,  AND  RECOMMENDATIONS 

Results 

A  new  and  effective  means  of  analyzing  aircraft  motion  has  been 
introduced  to  the  USAF  Test  Pilot  School.  These  methods  have  been 
used  by  selected  members  of  the  TPS  during  the  Have  Delay  Test  Project 
(Ref  10).  Aircraft  motion  and  aircraft  flight  control  systems  can  now 
be  examined  at  the  Test  Pilot  School  by  using  frequency  response 
techniques  using  a  program  called  FRAN.  In  addition,  complex  aircraft 
systems  which  have  been  analyzed  by  these  techniques  can  now  be 
reduced  to  Lower  Order  Equivalent  Systems  using  a  program  called  LOES. 
The  use  of  frequency  response  methods  for  pilot-in-the-loop  analysis 
and  the  use  of  lower  order  equivalent  systems  provide  the  USAF  Test 
Pilot  School  with  powerful  analysis  tools. 

The  main  result  of  this  thesis  effort  is  the  Test  Pilot  School 
student  now  has  the  resources  available  to  analyze  actual  aircraft 
response  in  the  frequency  domain.  Time  histories  of  control  inputs 
and  aircraft  responses  can  be  processed  to  produce  Bode  Plots  such  as 
shown  in  Figure  28.  Matching  of  actual  frequency  response  data  with 
Lower  Order  Equivalent  Systems  is  now  available  to  the  Test  Pilot 
School  for  a  larger  selection  of  Lower  Order  Transfer  Functions.  The 
Have  Delay  project  alone  utilized  the  software  generated  for  this 
thesis  for  27  different  configurations,  generating  27  seperate  sets  of 
magnitude  and  phase  plots,  power  spectral  densities  plots,  and  Lower 
Order  Equivalent  System  fits. 

Problems 

One  problem  associated  with  the  frequency  response  technique  is 
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that  the  method  is  not  designed  to  handle  non-linearities. 


Aerodynamic  non-linearities  are  kept  to  a  minimum  if  the  task  is 
designed  requiring  only  small  pertubations  in  angle  of  attack,  load 
factor,  etc.  Non-linearities  would  occur  if  the  pilot  re-trimmed  the 
aircraft  during  the  maneuver  thus  invalidating  the  stick  force  time 
history.  The  flight  control  system  can  be  a  large  source  of  non- 
linearities.  For  example  if  a  stick  controller  contains  a  breakpoint 
with  significantly  different  gradients  on  each  side  of  the  breakpoint, 
non-linear  effects  will  corrupt  the  analysis.  Parabolic  stick  force 
shaping  would  also  add  a  non-linear  contribution.  Non-linear  effects 
in  flight  controls  are  not  uncommon  and  frequency  analysis  methods 
which  attempt  to  account  for  these  non-linear  effects  need  to  be 
examined. 

It  appears  that  the  majority  of  the  uncertainty  in  the  frequency 
response  plots  generated  by  program  FRAN  can  be  attributed  to  the 
complete  lack  of  signal  processing  or  filtering  prior  to  the  signal 
being  digitized.  Gardenshirt,  Williams,  and  others  (Ref  15),  stress 
the  importance  of  properly  processing  and  filtering  flight  test  data 
prior  to  use.  Passing  the  data  through  a  reliable  low  pass  filter 
prior  to  digitizing  can  reduce  a  portion  of  the  aliasing  errors.  If 
the  original  waveform  was  not  sampled  frequently  enough,  then  the 
waveform  may  be  reconstructed  or  interpolated  after  digitizing  by  a 
number  of  reconstruction  techniques,  either  linear  or  higher  order 
interpolation.  Interpolation  can  also  introduce  a  number  of  errors  of 
m  the  original  waveform.  Methods  of  using  higher  sample 
upon  reconstructed  waveforms  and  use  of  anti-aliasing 


filters  are  not  in  use  by  the  USAF  Test  Pilot  School.  These 
techniques  need  to  be  examined  as  a  means  of  improving  results  of 
frequency  response  programs  (FRAN). 

The  manner  in  which  the  lower  order  equivalent  system  fit  is 
applied  to  the  frequency  response  data  has  been  shown  to  be  an 
accurate  and  reliable  manner  of  determining  equivalent  system 
parameters.  However  improvements  are  still  warranted.  The 
Proposed  Mil  Handbook  (Ref  8)  suggests  that  one  may  wish  to  use  a 
slightly  different  cost  function  and  may  even  want  to  use  different 
transfer  functions  for  the  LOES  fits.  The  industry  needs  to 
standardize  the  cost  functions,  frequency  range,  and  transfer 
functions  to  be  used  for  frequency  response  methods. 

While  much  improvement  in  noise  reduction  was  achieved  with  the 
use  of  averaging  overlapping  windows  and  the  use  of  faster  sampling 
rates,  the  frequency  response  produced  using  the  programs  developed 
for  the  Test  Pilot  School  are  still  too  noisy  for  analysis  at  the 
higher  frequencies  (10-20  rad/sec).  Analysis  of  flight  controls, 
especially,  require  accurate  information  past  10  radians  per  second 
which  is  not  currently  available.  Additional  studies  should  be  made 
into  ways  of  decreasing  noise  and  aliasing  effects  and  thus  increasing 
the  reliablity  of  the  Test  Pilot  School  new  frequency  analysis 
program. 
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FORTRAN  SOURCE  CODE  LISTING 
PROGRAM  FRAN 
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PROGRAM  FRAN 

C  xxxxxxxxxxxtxxxxxxxxxxxxxxxzxxxxxxxxxxxxxxxxxxx 
.  C  THIS  PROGRAM  WAS  DEVELOPED  TO  FULLFILL  THE 
C  REQUIREMENTS  OF  A  MASTERS  THESIS  BY  ALLAN  T.  REED 
C  DEVELOPED  AT  THE  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 
C  FROM  1984-8S  PROGRAM  USES  UP  TO  512  DATA  POINTS 
C  FROM  TWO  SIMULTANEOUS  PROCESSES  NO  DATA  PROCESSING 
C  (I.E.  NOISE  REDUCTION  OR  FILTERING)  IS  PERFORMED 
C  BY  THIS  PROGRAM.  DATA  PROCESSING  MUST  BE  DONE  PRIOR 
C  PROGRAM  DIVIDES  RAW  DATA  INTO  SEPARATE  ’WINDOWS*  FOR 
C  PROCESSING  THEN  PERFORMS  FFT  ON  THE  WINDOWED  DATA 
C  FREQUENCY  INFORMATION  IS  GENERATED  FOR  EACH  WINDOW 
C  AND  THEN  ALL  OF  THE  WINDOWED  DATA  IS  AVERAGED 
C  TOGETHER  TO  MININIMI2E  THE  EFFECTS  OF  NOISE  OR  ALAISING 
C  PROGRAM  OPTIMIZED  BY  LONG  TIME  HISTORIES  AND  FAST 
C  SAMPLE  RATES  (20  SAMPLES/SEC  MINIMUM) 

C  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
C  A&B  ARE  INPUT  COEFFICIENTS  C&D  ARE  OUPUT  COEFFICIENTS 
DIMENSION  AC 512) .BC512) . D(512) . RC512) . GXCS12) ,W(S12) 

DIMENSION  GYC512).  GXYCS12) .CC512) .SC512) 

DIMENSION  AT( 1024). CTC 1024) ,PA(S12).QA(512) 

BYTE  FILEC12) 

DATA  FILE/8X0. * . * . 'D' . ’A* . 'V / 

WRITE  (S.100) 

100  FORMAT  C’  ENTER  NUMBER  OF  DATA  POINTS  AND  POWER  OF  TWO  ’ ) 
READC5. 10DN.NU 

101  FORMAT  (14) 

PI -3. 1415926535 

c  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxsx 
C  INITIALIZATION  FOR  THE  WINDOW  SIZE 
ND2-NX2 
ND4-N/4 
ND8-N/9 
NUM2-NU-2 
C 

WRITE  (5.270O) 

2700  FORMAT ( ’  ENTER  DATA  SAMPLE  RATE  IN  SECONDS.  I.E.  .05’) 

READ  (5.102)  SR 
C 

DO  1  L-1.ND2 

GX(L)-0.0 

GXY(L)-0.0 

PA(L) “0. 0 

QA(L)-0.0 

GY(L)-0.0 

C  W(L)  IS  THE  ARRAY  OF  FREQUENCIES 
W(L)-(LX2.0*PI )/( FLOAT (N)XSR) 

W(L)-ALOG10(W(L)) 

1  CONTINUE 

C  PROGRAM  USES  A  MINIMUM  50*  OVERLAPPING  WINDOW  FOR  FFT  CALCULATIONS 
WMIN-0.01 
WMAX-100. 

WRITE(5>  200) 

200  FORMAT ( ’  ENTER  NAME  OF  INPUT  DATAi  NNNNN.DAT  ’) 

ACCEPT  52.  (FILEC I). 1*1.0) 

52  FORMAT (8A1) 

OPEN(UNIT-l. NAME-FILE. TYPE- ’UNKNOWN’ ) 

WRITEC5.40) 

40  FORMAT ( ’  READING  IN  INPUT  DATA  ’,/) 

DO  41  I-l.N 

C  ARRAY  ATCI)  IS  A  FILE  FOR  ENTIRE  INPUT  DATA  ARRAY 
READC1. 102)AT( I ) 

41  CONTINUE 
WSITEC5. 201 ) 
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201  FORMAT ( *  ENTER  NAME  OF  RESPONSE  DATA'  t'NHN.  DAT  *  ) 

ACCEPT  52. (FILEC I). 1-1.8) 

CLOSE  (UNIT-1) 

OPEN  (UNIT-2. NAME-FILE. TYPE- ’ UNKNOUN’ ) 

WRITECS.241) 

241  FORMAT ( ’  READING  RESPONSE  DATA*  ,/) 

DO  42  I-l.N 

C  ARRAY  CT( I )  IS  A  FILE  FOR  ENTIRE  OUTPUT  DATA  ARRAY 
R£AD(2. 102)CT ( I ) 

42  CONTINUE 

C  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

c 

I  F(N.  EQ .  512)  KWNDOW-7 
IF(N.EQ. 1024)  KWNDOW-15 

C  A  TOTAL  OF  7  OR  15  WINDOWS  WILL  BE  AUERAGED  (50*4  OUERLAP) 

DO  2  K-l.KWNDOU 
URITEC5.50) 

50  FORMAT ( '  WINDOWING  DATA  PLUS  BLACKMAN  FILTERING  '  ./) 

C  ASSIGN  INPUT  AND  OUTPUT  DATA  FOR  EACH  FFT  WINDOW 
C  J  -  START  POINT  OF  CURRENT  WINDOW 

C  J1  •  STOP  POINT  OF  CURRENT  WINDOW 

C 

J*(CK-1)*64)  +  1 

Jl-J+128-1 

M-0 

C  LOAD  APPROPRIATE  AMOUNT  OF  POINTS  AND  BLACKMAN  FILTER 
C  BLACKMAN  FILTERING  IS  SIMILAR  TO  COSINE  OR  HANN  FILTERING 
C 

DO  3  I-J.J1 
M-M+l 

R(M) *  AT( I ) 

S(M)«  CT( I ) 

ARG-FLOAT(M)/FLOAT( 128) 

R(M)-R(M)*C.42-.50*COSC2.*PI*ARG)+.08*COS(4.*PI*ARG)) 
S(M)-S(M)*(. 42- .  50*COS  ( 2 . *PI*ARG)+. 08*COS(4. *PI*ARG) ) 

3  CONTINUE 

102  FORMAT  CF15.9) 

WRITEC5.53) 

53  FORMAT (’  COMPUTING  FFT  ON  INPUT  AND  OUTPUT  DATA’./) 

CALL  FFTCRfS. 128.7. 1) 

C 

C  FFT  CONUERSION  COMPLETE 
C 

c  XXXXXXXXXXXXXXXXXX XXXXX X X XXX xxxxxxxxxx 

c 

c  SEPERATE  FFT  DATA  INTO  REAL  AND  IMAGINARY  COMPONENTS  OF 
C  INPUT  AND  OUTPUT  SIGNALS. 

C 

DO  4  1-1,64 
L-I+l 

NMI-128-I+1 

AC  I )-(R(L)+R(NMI ))/2. 

B(  I  )  ■  ( S ( L ) -S ( NM I ) ) 42 . 

C( I )-(S(L)+S(NMI ) )/2. 

DC  I  )—  (R(L)-R(NMI  ))/2. 

4  CONTINUE 
C 

WRITEC5. 104 )J . J1 

104  FORMAT ( ’  COMPUTING  PARAMETERS  FOR  WINDOW  UALUES* . 14, ■  TO’, 14. 4 
DC  5  1-1.64 

P-A( I )*C( I )  ♦  B( I )XD( I ) 

Q-A( I )*D( I )  -  B( I  1*CC I ) 

C  COMPUTE  CROSS  SPECTRAL  DENSITY  OF  INPUT  TO  OUTPUT.  SXY 
GXYCD-  GXY (  I )  +SGRT ( P**2 .  +  0**2.) 

C  COMPUTE  POWER  SPECTRAL  DENSITY  OF  INPUT.  GX 
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GXC  I)-GXC  I)+A(  I  )XX2.  +  BCI)*X2. 

pac  i )-pac  n  +  p 

QAC I )-QAC I  )  +  Q 

C  COMPUTE  GY  •  PSD  OF  OUTPUT  SIGNAL,  GY 
GYC  I  )-GYC  I  )-HDC  I  )XX2.  +DCI)*X2. 

5  CONTINUE 

2  CONTINUE 

C 

C  tXtXXXXTttXXK*X*XKXXXXXt*XXX%XXXXXXX*X% 

C 

CLOSE (UN  IT-2) 

C  AUERAGE  ALL  UALUES 
UR  I TE  ( 5 , 60 ) 

60  FORMAT ( ’  AUERAGING  PARAMETERS  AND  INTERPOLATING  FREQUENCIES'./) 

MID  -64 
LAST- 129 
NTL-LAST-1 
F  -  FLOAT  ( K14NDOU ) 

DO  6  I-l.MID 
GXYC I )-GXYC I )/F 
GXC I )-GX( I )/F 
GYC I )-GY( I )/F 
PA( I ) -PA( I )/F 
QA( I )-QAC I )/F 

6  CONTINUE 

C  ASSIGN  TEMPORARY  UALUES  TO  FILL  OUT  FREQUENCY  RANGE 
90  DO  66  I-l.MID 

AC  I ) -GXYC I ) 

B( I )-GX< I) 

C( I ) -GYC I ) 

ATC  D-PACI) 

cTcn-QAcn 
66  CONTINUE 

C  FILL  OUT  ENTIRE  FREQUENCY  RANGE 
DO  7  1-2. LAST, 2 
ID2-I/2 
GXYC 1 3 -AC ID2) 

GXCD-BCID2) 

GYC I )-CC ID2) 

PAC I ) -ATC ID2) 

QAC I ) -CT( ID2) 

7  CONTINUE 

C  INTERPOLATE  INTERMEDIATE  UALUES 
DO  9  I-3.NTL.2 
IM1-I-1 
IP1-I+1 

GXYC  I ) -(GXYC  IM1)+GXYC  IF1)  )/2. 

GXC  I  )-(GXC  IMD+GXC  I  PI)  )/2. 

GYC  I  )  •  C GY (  IMD+GYC IP1)  )/2. 

PAC  D-CPAC  IMD+PAC  I  PI)  )/2. 

QAC  I)-CQA(  IMD+QAC  IP1)  )/2. 

0  CONTINUE 

C  LINEARLY  INTERPOLATE  FIRST  AND  LAST  UALUES  OF  RECORDS 

GXY ( 1 ) • C  GXY  C  2 ) -GXY  C  3 ) ) +GXY ( 2 ) 

GXC  1)-CGX(2)-GX(3)  )-H3XC2) 

PAC 1 ) - C  PAC2)-PA(3 ) ) +PAC2) 

QA  C 1 ) - C  QA  C  2 ) -QA  C  3 ) ) +QA ( 2 ) 

GY ( 1 )  *  ( GY ( 2 ) -GY C 3 )  )-K5YC2) 

NN-NTL-1 

GXY  ( LAST )  •  C  GXY  C  NTL )  -GXY  C  NN )  )  +GXY  ( NTL ) 

GX ( LAST ) - C  GX  C  NTL ) -GX ( NN ) 5+GXCNTL) 

PA C LAST) • C PA C NTL ) -PA C NN ) ) +PA C NTL ) 

QA C LAST )• C QA C NTL ) -QA C NN ) ) +QA C NTL ) 

GY C LAST ) • ( GY ( NTL ) -GY C NN ) ) -H3Y C NTL ) 

C  IF  ND2  DATA  POINTS.  FINISHED!  IF  NOT.  REPEAT  DOUBLING  ROUTINE 
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IF  C  MID  .  EQ.  ND4)  GO  TO  99 
LAST  -  LAST *2 
NTL  •  LAST-1 
MID  -  M1DX2 
GO  TO  98 
99  CONTINUE 

C  COMPUTE  HXY-GXYxGX  AND  GAMMA  SQUARED  UALUES 

OPEN (UNIT-1. NAME-’ BODE. BDE’ .TYPE-’ UNKNOWN’ ) 

C  COMJERT  HXY  TO  DECIBLES 
C  LET  A(I)  BE  THE  HXYC I)  UALUES 

C  LET  B( I)  BE  THE  GAMMA  SQUARED.  GAMMA2C I )  UALUES 
C  LET  CCI1  BE  THE  PHASE C I )  UALUES 
DO  9  I-1.ND2 
AC  D-GXYC I )/GX( I ) 

AC  I )-20.*ALOG10(A( I ) ) 

BC I ) «20. X(GXY( I )*X2. )/(GX( I )*GYC I ) ) 

CCI) -ATAN2CQA( I ) . PA( I) ) *57. 29570 
IFCCCD  .GT.  90.)  C(I)-C(i)-360. 

C  COMJERT  GX  AND  GY  TO  POWER  SPECTRAL  DENSITIES 
GXC I )-10.*ALOG10(GXC I) ) 

GYC I ) « 10. *ALOG10( GY ( I ) > 

9  CONTINUE 

0PEN(UNIT-2.NAME-'GAMMA2.BDE’ .TYPE -’UNKNOWN’ ) 
WRITEC2, 1 3600 )ND2.UM IN. UMAX 
13600  FORMATCX. I3.G16.6.G16.6) 

DO  10  I-1.ND2 

URITEC2. 13700 )B( I ) ,C( I ) » W( I ) 

13700  FORMAT (X.3G16. 6) 

10  CONTINUE 
CLOSE  (UNIT-2) 

C  URITE  REMAINING  DATA  TO  DISK 

OPEN (UN IT-2. NAME-' GX.BDE’ .TYPE-’ UNKNCMT  ) 
WRITE(2. 13600)ND2. UMIN. UMAX 
URITE(1. 13600)ND2.LMIN.LT1AX 
DO  13  I-1.ND2 

WRITEd.  13700)  A(  I ) .  C(I).  U(I) 

WRITE(2. 13700)  GX(I).  C(I).  W(I) 

13  CONTINUE 
STOP 
END 
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APPENDIX  B 


FORTRAN  SOURCE  CODE  LISTING 
PROGRAM  LOES 
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c 

c 


PROGRAM  LCES  FOR  SYSTEM  IDENTIFICATION  BY  MATCHING  TRANSFER 
FUNCTION  TO  GIVEN  SET  OF  FREQUENCY  DATA.  TRANSFER  FUNCTION 
IS  SELECTED  BY  USER  THEN  PROGRAM  PERFORMS  CONJUGATE  GRADIENT 
C  SEARCH  TO  OPTIMIZE  COST  FUNCTION.  PROGRAM  REQUIRES  TWO  SUB- 
C  ROUTINES'  CFG  AND  CGFM  THESE  SUBROUTINES  PERFORM  THE 
C  COST  FUNCTION  CALCULATIONS  AND  THE  GRADIENT  SEARCH. 

PROGRAM  LOES 

DIMENSION  BCIS).G(IS). I  DC  IS) . AMP (50) .PHAC50) .  U(S1) .H(30) 
DIMENSION  URCS0? .DBC50) 

COMPLEX  PLANT C 50) 

EXTERNAL  CFG 
LOGICAL  KG 

C  .TRUE.  KG  MEANS  TO  COMPUTE  GRADIENT 

COMMON/ I DTN/PLANT.U. ID.K1 . INST. ND. KNT. NP. UR 
COMMON/ FMC/KOUNT . KG 

C  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
C  PROGRAM  VAR I  ABLE  LIST 
C 

C  «NP  SYSTEM  FREQUENCY  RESPONSE  MAGNITUDE  (UECTOR) 

C  B  MODEL  PARAMETERS  C INPUT  UAR I  ABLE  UECTOR) 

C  CFG  SUBROUTINE  UHICH  COMPUTES  COST  FUNCTION  AND  GRADIENT 

C  EPS  PARAMETER  CHANGE  DEFINING  SEARCH  CONUERGENCE  10E-5 

C  F  COST  FUNCTION  F-JCB) 

C  G  COST  FUNCTION  SCALED  GRADIENT  (UECTOR) 

C  H  STORAGE  UECTOR 

C  ID  INTEGER  UECTOR  THAT  IDENTIFIES  CORRESPONDING  B  PARAMETER 
C  IER  SEARCH  CONUERGENCE  PARAMETER  (0-CONUERGED '  SIZE<EPS) 

C  INST  BRANCHING  INSTRUCTIONS 

C  ITER  SET  EQUAL  TO  LIMITC200)  UNLESS  INST  •  4  OR  S 
C  J  INDEX  OF  ELEMENT  IN  UECTOR 

C  JD  NUMBER  OF  NONITERATIUE  CONSTANTS  IN  MODEL 
'  C  KG  LOGICAL  UAR I ABLE  .TRUE.  MEANS  COMPUTE  GRADIENT 

C  KNT  COUNT  OF  COST  FUNCTION  EVALUATIONS 
C  KOUNT  COUNT  OF  LINE  SEARCH  ITERATIONS 

C  K1  EXPONENT  OF  FREE  S  IN  MODEL  (INPUT  UAR I ABLE) 

C  LIMIT  MAXIMUM  NUMBER  OF  ITERATIONS  (SET  -  200) 

C  N  NUMBER  OF  ITERATIUE  PARAMETER  OF  B  UECTOR 

C  ND  NUMBER  OF  FREQUENCY  POINTS  OVER  WHICH  INTEGRATION  IS  DONE 

C  NP  TOTAL  NUMBER  ITERATIUE  PLUS’  CONSTANT  MODEL  PARAMETERS 
C  PHA  SYSTEM  FREQUENCY  RESPONSE  PHASE  (DEG)  (INPUT  UECTOR) 

C  PLANT  CONVERTED  SYSTEM  RESPONSE  TO  EQUIVALENT  COMPLEX  NUMBER 

C  SIZE  PARAMETER  STEP  SIZE  (SET  TO  0.1  AT  START  OF  SEARCH) 

C  THETA  SCRATCH  VARIABLE  FOR  RADIAN  PHASE 
C  W  RADIAN  FREQUENCY  VECTOR  (RADIANS/SECOND) 

C 

C  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
PI-3. 141592653S 
10  WRITE  (5,110) 

110  FORMAT ( ’  ENTER  NUMBER  OF  FREQUENCY  DATA  POINTS  FOR  CURVE  FIT’./) 
READ  (S.lll)ND 

111  FORMAT  (15) 

WRITECS. 1125 

112  FORMAT ( ’  FREQUENCY  DATA?  1-HAND  ENTER  2-FROM  FILE  FREQ. DAT  ’.// 
READ ( 5 , 140)  KK 

IF(KK  .EQ.  2)  GO  TO  15 
WRITECS. 103) 

103  FORMAT (*  ENTER  ALL  FREQUENCIES  IN  ORDER  (RAD/SEC) 

DO  2  I-l.ND 


READ  (5. 113)  U(  I ) 

2  CONTINUE 

WRITE  CS. 114) 

114  FORMAT C  ENTER  CORRESPONDING  AMPLITUDES  (DECIBLES)  ’,//) 


AD-A168  £04  EVALUATION  OF  A  FREQUENCY  RESPONSE  TECHNIQUE  FOR  2/2 

AIRCRAFT  SYSTEM  IDENTIFICATIONS)  AIR  FORCE  INST  OF 
TECH  HRIGHT-PATTERSON  AFB  OH  SCHOOL  OF  ENGI  A  T  REED 
UNCLASSIFIED  31  OCT  85  AFIT/GAE/AA/85J-2  F/G  1727  NL 


C  WRITE  (5. 103) 

193  FORMAT  (’  ENTER  MAX  NUMBER  OF  I TER I AT IONS  DESIRED’.//) 

C  READ  (S.lll)  LIMIT 

C  WRITE  (5.184) 

184  FORMAT  (’  ENTER  UALUE  OF  EXPONENT  OF  THE  FREE  S  *.//) 

C  READ  (5. Ill )  K1 

C  WRITE  CS.18S) 

185  FORMAT  (’  PROGRAM  WILL  ITERATE  ON  ALL  B  UALUES  UNLESS  YOU  ’./. 

1*  CHOOSE  THE  LAST  ENTERED  B  UALUES  TO  BE  FIXED.  ’./. 

2*  ENTER  NUMBER  OF  FIXED  B  UALUES  (NORMALLY  -  0)  ’.//) 

C  READ  (S.lll)  JD 

C  xxxxxxxxxxtxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

36  N-N-JD 

NP-N+JD 

C  NP  IS  THE  TOTAL  NUMBER  OF  TERMS  BOTH  ITERATIUE  AND  CONSTANTS 
WRITE  (5.130)  LIMIT. N. Kl. JD. ( ID(J) . J-l.NP) 

130  FORMAT ( ’  LIMIT. N.Kl. JD-’ .413. *  ID-’15I3) 

GO  TO  87 

777  WRITE  (5.186) 

186  FORMAT(/. ’  INSTRUCTIONS?  1-  SEARCH  2-  NEW  FREQUENCY  DATA  ’./. 

1'  3-  NEW  TRANSFER  FUNCTION  4-  PRINT  &  SAUE  RESULTS  5-STOP  ’) 

READ  (5.140)  INST 
140  FORMAT  (ID 

GO  TO  (40,10.30,50.1000). INST 
50  URITE  (5.160) 

160  FORMAT ( 6H  MODEL, /. 8X. 1HW. 14X.3HAMP. 11X.3HPHA) 

40  SIZE*. 1 

C  SIZE  IS  THE  PARAMETER  STEP  SIZE 
EPS-1. E-6  _ 

C  EPS  IS  THE  PARAMETER  CHANGE  DEFINING  CONUERGENCE 
I TER-LIMIT*(1- INST/4) 

C  IF  INST  -  4  ITER  BECOMES  ZERO  AND  NO  ITERATIONS  TAKE  PLACE 
KNT-0 

C  KNT  RECORDS  NUMBER  OF  TIMES  THE  COST  FUNCTION  IS  EUALUATED 
CALL  CGFMtCFG.N.B.F. G.SIZE.EPS. ITER. IER.H) 

WRITE (5. 150)  F. IER.KOUNT, KNT. SIZE 
C  F  IS  THE  COST  FUNCTION 

150  FORMATC  F-’F15.5./. 

1*  ERROR  MSG  -  (IER)  (0-CONUERGE. 1-NO  CONVERGE, 2 -ERROR)  ’ ,/. 

2’  IER-’. 12.’  KOUNT. KNT. SIZE-’ .13.14. 1PE10.3) 

87  CONTINUE 

WRITE(S, 151)  (B(J) . J-l.NP) 

151  FORMATC  B-’.8F11.4) 

GO  TO  777 

1000  STOP 
END 


SUBROUTINE  CFG  COMPUTES  THE  COST  FUNCTION  AND  GRADIENT 
SUBROUTINE  CFG(N.B.F.G) 


CFG  COMPUTES  THE  COST  FUNCTION  AND  GRADIENT 
DIMENSION  B( 1) ,G( 1 ) • !D(15) , W(S1 ) 

COMPLEX  PLANTC50).Z(15).S.ST.GM.E 
LOGICAL  KG 

COMMON/ 1 DTN/PLANT.W.  ID.K1.  INST. ND. KNT. NP 
COFMON/FMC/KOUNT . KG 

xxxxxxxxxxxxxxxxxxxxxxxxxxxtxxxxxxxxxxxxxxxxxxxx 

SUBROUTINE  CFG  UARIABLE  LIST 

AM  MAGNITUDE  OF  G(JW.B) 

DEG  ANGLE  OF  G(JW,B)  IN  DEGREES 
E  ERROR  BETWEEN  MODEL  AND  DATA 

GM  PARTIAL  PRODUCT  BECOMES  G(S.B) 


mi 


r  -  v  -  -  '  ►  »  >  - 


UALUE  OF  IDUO 
FREQUENCY  INDEX 
PARAMETER  INDEX 
S 

TEMPORARY  UALUE 
SAUED  W  UALUE 

PARTIAL  PRODUCT  IN  MODEL  GRADIENT  (UECTOR) 


****************************** ******************* 

I F C I NST  .EQ.  4)  0PEN(UNIT-2.NAME«’RESULTS. DAT’ .TYPE-' 
KNT-KNT+1 
DO  30  J-l.N 
G( J3 -0. 

CONTINUE 

ws-uc i 3 

F-0. 

DO  100  J-l.ND 
S-CMPLXC0. ,U(J)) 

S  -  LAPLACE  UARIABLE 
GM-S**K1 
DO  0  K-l.NP 
IDK-IDCK) 

GO  TO  (1.3. 1. 2.8.2. 8.4). I DK 
START  HERE  IF  DEALING  WITH  A  ZERO  OR  A  POLE 
GM-GM*(S/B(K)+1.  )*X(2-1DK) 

GM  IS  A  PARTIAL  PRODUCT  OF  TRANSFER  FUNCTION 
IF  (KG)  Z(K3-S/(S+S(K33*(FLOAT(IDK)-2. ) 

GO  TO  0 

START  HERE  IF  USING  MODEL  AS  PLANT 

ST-(S/B(K+1)3**2+2.*S*B(K)/BCK+1)+1. 

ST  IS  A  TEMPORARY  UALUE 
GM-GM*ST**(S-IDK)  • 

IF(.NOT.KG)  GO  TO  0 

Z(K)-2.*S*B(K)/<B<K+1)*ST)*(5.-FL0AT(IDK)> 

Z<K+1)— Z(K)*1.+S'((B(K)*B(K+1))) 

GO  TO  8 

Z  IS  THE  PARTIAL  PRODUCT  IN  MODEL  GRADIENT 
START  HERE  IF  DEALING  UITH  THE  GAIN 
GM-GM*B(K) 

Z(K)-CMPLX<1. ,0. 3 

GO  TO  8  _ 

START  HERE  IF  DEALING  WITH  TRANSPORT  LAG  (DEAD  TIME) 
GM-GM*CEXF(-S*BCK) ) 

ZCK)— SXBCK) 

CONTINUE 

E-GM-PLANT(J) 

E  REPRESENTS  THE  ERROR  BETWEEN  THE  MODEL  AND  THE  PLANT 
F*F+(W(J+1 )-US)*R£AL(E*CONJG(E) * 

IF( .NOT.  KG)  GO  TO  100 
ST • ( U ( J+ 1 ) -WS ) *E*CONJG ( GM ) 

DO  80  K-l.N 

G ( K ) - G ( K ) + REAL ( ST*CONJG ( Z ( K ) ) ) 

CONTINUE 

I F ( INST. EQ. 1 3  GO  TO  100 
AM-CABS(GM) 

DB-20.*ALOG10(AM) 

AM  IS  THE  MAGNITUDE  OF  THE  PARTIAL  PRODUCTS 
DEG - ATAN2 ( A 1 MAG ( GM 3 , REAL ( GM ) 3 *S7 . 29S780 
IF  (DEG  .GT.  90. 3  DEG -DEG- 360. 

WRITE(2» 120)  U( J) . DB. DEG 
WRITE(S. 1203  U( J 3 , DB. DEG 
0  FORMAT (3F15. 5) 

IF( INST. EQ. S3  PLANTc J3-GM 
0  WS-W(J) 

I F ( I NST  .EQ.  43  URITE(2.2003F. (B(J).J-l.NP) 


UNKNOWN’ 3 


'vvVv' v'v^v'* v 


200  FORMATt’  F-* .F15.5./. ’  B-’.9F11.4) 

IFCINST  .EQ.  4)  CLOSE (UN IT-2) 

RETURN 

END 

C  txxxxxxxxiixxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx. 

C  SUBROUTINE  CGFM  CONDUCTS  THE  CONJUGATE  GRADIENT  SEARCH 
C 

SUBROUTINE  CGFM(CFG. N.B.F. G.SI2E.EPS. ITER. IER.H) 

C 

C  CGFM  PERFORMS  THE  ACTUAL  MINIMIZATION  OF  THE  CONJUGATE 
C  GRADIENT  SEARCH  AS  STARTED  BY  SUBROUTINE  CFG 
DIMENSION  B(1).H(1) ,G(1) 

LOGICAL  KG 
COMMON/ FMC/KOUNT . KG 

C  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
C  SUBROUTINE  CGFM  UARIABLE  LIST 

c 

C  BETA  CONJUGATE  DIRECTION  WEIGHTING 

C  FS  SAUED  F 

C  FSS  SAUED  FS 

C  J  PARAMETER  INDEX 

C  JN  J/N 

C  K  INDICATOR  FOR  STEP  SIZE  REDUCTION 

C  L  NUMBER  STEP  SIZE  INCREASES  WITHIN  ITERATION 

C  NCYC  NUMBER  OF  ITERATION  BEFORE  RESTARTING  CONJUGATE  SEARCH 

C  R1.R2.R3  TERMS  IN  QUADRATIC  CURUE  FIT 

C  SCALE  STEP  SIZE  SCALE  FACTOR 

C  SCASU  SAUED  SCALE 

C  STEP  STEP  SIZE 

C  TSAUE  SAUED  TSQR 

C  TSQR  SQUARED  GRADIENT  TERMS  SUM 

C  XI. X2  PARTIAL  PRODUCT 

C 

C  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
C  INITIALIZATIONS 
KQUNT-0 
STEP-2. 

IER--5 
5  BETfl-0 

C  BETA  IS  THE  WEIGHTING  ON  THE  CONJUGATE  DIRECTION 
NCYC-0 

C  NCYC  IS  THE  NUMBER  OF  ITERATIONS  BEFORE  RESTARTING  CONJUGATE  SEARCH 
IS  K-0 

DO  20  J-l.N 
H(J)-B(J) 

20  CONTINUE 

KG-.TRUE. 

CALL  CFGCN.B.F.G) 

C  TEST  FOR  STOPPING  SEARCH 

IFOCOUNT  .GE.  ITER)  IER-1 
IFCSIZE  .LT.  EPS)  IER-0 
40  I F ( IER  .GT.  -2)  RETURN 

C  COMPUTE  PAST  GRADIENT  WEIGHTING 
TSQR-0. 

C  TSQR  IS  THE  SUM  OF  THE  SQUARED  GRADIENT  TERMS 
DO  50  J-l.N 

50  TSQR-TSQR-H3(J)**2 

IFCNCYC  .EQ.  0)  GO  TO  60 
BETA  -  TSQR/TSAUE 
60  SCALE -0. 

DO  70  J-l.N 
JN-J+N 

H( JN)--G( J)+BETA*H( JN) 

70  SCALE -SCALE+ABS ( H ( JN ) ) 

I F( SCALE  .GT.  0. )  GO  TO  90 


87 


00 


GO  TO  40 
SCALE-SIZE/SCALE 
TSAUE-TSGR 
NCYC-NCYOl 
FSS-F 
L-l 

C  L  COUNTS  NUMBER  OF  ITERATIONS 
C  UPDATE  THE  UALUES  OF  THE  B’S 
100  DO  110  J-l.N 
JN-J+N 

110  B(J)-H(J)*ABSU.+SCAL£*HCJN)) 

FS-F 

KG-. FALSE. 

CALL  CFG(N.B.F.G) 

C  LOGIC  CHANGES  LINE  SEARCH  STEP  SIZE  OF  CONCLUDES  SEARCH 

IFCK. GT.0)  GO  TO  120 
IF(F.LT.FS)  GO  TO  130 

IFCL. GT. 1)  GO  TO  140 
GO  TO  1S0 

120  IF(F.LT.FSS)  GO  TO  140 
GO  TO  160 

130  SCASU  -  SCALE 

C  SCASU  IS  THE  SAUED  UALUE  OF  THE  SCALE 
L«L>1 

SCALE -SC ALEXSTEP 
FSS-FS 

IFCL.LT.lS)  GO  TO  100 
IER  -2 

C  IER  UALUE  OF  2  MEANS  ERROR  HAS  OCCURRED 
GO  TO  40 

C  NOW  FIT  QUADRATIC  CURUE  TO  3  PTS.  BRACKETING  LINE  SEARCH  MIN 
140  DO  140  J-l.N 
JN-J+N 
Rl-H(J) 

R2-H(J)*ABS<1.+SCASU*H<JN)) 

R3-BC J)  _ 

IF(L.GT.S)  R1-R1+<R2-R1)/STEP 
X1«(FSS-FS)*(R1-R3) 

X2- ( FSS-F )*(R2-R1) 

IF(ABSCR2-R3).GT.EPS^4.)  SO  TO  147 
IF(L.GT.l)  B( J)-R2 
GO  TO  140 

147  BCJ)-(X1*(R1+R3)4X2*(R14R2))/'((X1+X2)*2. ) 

140  IF(B(J)*H(J).LE.0. )  B( J)— 1 . *B( J)+EPS*R3 

C  UPDATE  SEARCH  VARIABLES 

SIZE-SIZE*CFL0ATCL)+2.  )✓ 4. 

KOUNT-KOUNT+1 
IFCNCYC.GT.N)  GO  TO  5 
GO  TO  15 

150  SCASU-SCALE 

SCALE-SCALE/ ( 1 . +STEP ) 

K-K+l  _ 

SIZE-SIZE/( 1 . +5TEP) 

GO  TO  100 

160  SIZE-SI  ZLs ( 1 . +STEP ) 

DO  100  J-l.N 
100  B( J)-H(J) 

GO  TO  5 


PROGRAM  F410 

F4  AT  FLIGHT  CONDITION  NO.  i0 

DIMENSION  A(4.4).  CMATC4, 4) . EAT(4. 4) . EATM1(4. 4) 

DIMENSION  IM(4.4).AT(4.4). ATT2 (4.4). FCAPT (4.4) 

DIMENSION  BU(4.4).X(4.4),AI(4,4),ASQ(4,4),DEL(512) 

DOUBLE  PRECISION  CHAT 

REAL  IM 

INTEGER  RA.CA 

BYTE  FILE(12).0FILE(14) 

DATA  FILE/0X0. * . * . ’D’ . ’A* , ’T'/ 

RA-4 

CA-4 

0  FORMAT  (ID 
NUT-0 

INITIALIZE  ALL  MATRICIES 
DO  1  I  -  l.RA 
DO  2  J  -  l.RA 
ASQl I.J)-0.0 
EAT( I . J) -0.0 
EATM1( I. J)-0.0 
AI( I,  J)-0.0 
ATC I . J)-0.0 
X(  I.J)-0.0 
ATT2( I. J)-0.0 
FCAPT ( I • J)-0.0 
BU(I,J)-0.0 
IM(I.J)  -  0.0 
CONTINUE 

SYSTEM  A-MATRIX  HERE 


A(  1  ■ 

1) 

■— 

.00520 

A(l. 

2) 

•— 

.00501 

A(  1 

3) 

•  — 

100.3 

A(l. 

4) 

■  — 

31.947 

A(2. 

1) 

■  m 

000474 

A(2. 

2) 

■— 

.319 

A(2, 

3) 

-1730 

A(2. 

4) 

-1 

.842 

A(3. 

1) 

00175 

A(3. 

2) 

.0133 

A(3, 

3) 

■  — 

.3130 

A(3, 

4) 

-0 

.0 

A(4, 

1) 

-0 

.0 

A(  4, 

2) 

-0 

.0 

A(4, 

3) 

-1 

.0 

A(4, 

4) 

-0 

.0 

WRITE  (5.13) 

FORMAT  ("ENTER  TIME  OF  OUTPUT .  STEP  SIZE. AND  PRINT  STEPS 
READ  (5.200)  TIME.CAPT 
READ  (5.201)  IPRNT 
FORMAT  (115) 

FORMAT  (1F15.9) 

STIME-0. 0 
DTI ME  ■  TIME+.l 
Bll-3.42 
B2 1—20.7 
B31— 4.90 
B41-0.0 
WRITE  (5.203) 

FORMAT  (’  ENTER  NUMBER  OF  DATA  POINTS  ',/) 

READ  (5.201)  K 
WRITE  (5.51) 

FORMAT ( ’  ENTER  NAME  OF  INPUT  DAT'  NNNNNN.DAT’ .✓) 

ACCEPT  52.  (FILE(t), I-l.B) 

FORMAT  (0A1) 


OPEN (UN IT-3 .NAME- FILE. TYPE -'UNKNOWN* ) 

DO  5  I-l.K 
READC3.200)  DEL( ! ) 

S  CONTINUE 

CLOSECUNIT-3) 

C  OPEN  FILES  FOR  OUTPUT 

OPEN ( UN I T-l. NAME-’ UELOCITY.DAT' .TYPE- ’UNKNOWN’ ) 

OPEN CUN IT-2. NONE -*OLPHO. DOT' .TYPE -’UNKNOWN* ) 

OPEN ( UN IT-3. NONE- *Q. DOT’ .TYPE- ’UNKNOWN* ) 

OPENCUNIT-4, NONE- ’THETO. DOT’ .TYPE- ’UNKNOWN’ ) 

C  START  BY  SETTING  01  EQUAL  TO  0  MATRIX 
DO  3  I-l.RO 
DO  4  J-l.Cfl 
4  01 (I , J)-0( I . J) 

3  CONTINUE 

C  1NUERT  THE  0  MATRIX  USING  INUSUBl  A  MATRIX  INUERSION  SUBROUTINE 
C 

CALL  INUSUB  (A.AI.RA.CA) 

C 

C  MOTMUL  IS  A  MATRIX  MULTIPLICATION  SUBROUTINE 
CALL  MOTMUL  C A, A. CMAT. RA) 

DO  31  2-1. RA 
DO  32  J-  l.CA 

32  OSQ  (I.J)  -CMAT  (I.J) 

31  CONTINUE 

JJ-0 

C  EAT  -  I  +  AT  +  (0X*2)TXX2  /  2  FACTORIAL  +  (AX*3)TX*3  x  3  FACTORIAL 
7  CONTINUE 

DO  91  I-l.RA 

DO  92  J-l.CA 

AT  (I.J)-  A(I,J)*CAPT 

ATT2(  I .  J)-  (ASQ(I.J)*CAPT**2. )/2. 

92  EAT (I.J)  -  ATd.J)  +ATT2(I.J) 

91  CONTINUE 

DO  64  I-l.CA 
EAT  (I.I)-EAT  (1. 1)4-1. 0 
64  CONTINUE 

CALL  MOTMUL  (A. ASQ.CMAT. RA) 

DO  33  I-l.RA 
DO  34  J-l.CA 

34  EAT(I,J)-EAT(I.J)4-(CMATCI.J)*CAPT**3.  )x6.0 

33  CONTINUE 

C  FIND  F(T)  MATRIX  FUNCTION 
DO  19  I-l.RA 
DO  29  J-l.CA 

29  EATMl(I.J)-EATU.J) 

19  CONTINUE 

DO  21  I-l.RA 

EATMl(I.I)  •  EAT(I.I)  -  1.0 

21  CONTINUE 

22  JJ-JJ+l 

BUC1, l)-BllXDEL(JJ) 

BU(2> 1)-B21*DEL(JJ) 

BU<3. 1 ) -B31XDELC JJ) 

CALL  MOTMUL  (EATM1.BU.CMAT.RA) 

DO  23  I  -  l.RA 
DO  24  J  •  l.CA 

24  FCAPTCI.J)  •  CMAT (I.J) 

23  CONTINUE 
C 

CALL  MOTMUL  (AI.  FCAPT . CMAT . RA ) 

DO  2S  I  •  l.RA 
DC  26  J-  l.CA 
26  FCAPT (I.J) -CMAT (I.J) 

25  CONTINUE 


91 


77  CONTINUE 

C  MAIN  ITERATION  OCCURS  HERE 

CALL  MATMUL  (EAT. X. CHAT. RA) 

DO  37  I-l.RA 
DO  38  J-l.CA 

38  X( I ,  J)  -  CMATC I . J) 

37  CONTINUE 

990  CONTINUE 

DO  39  I-l.RA 
DO  41  J-l.CA 

41  X(I.J)-XCI.J)  +FCAPT  C I » J ) 

39  CONTINUE 

IFCSTIME  .EQ.  0.0)  URITEC5.il) 

11  FORMATC//.'  DELTA  UALUES * . SX , * UELOC I TY ’ . 8X . ' ALPHA ’ , 8X 
l'P  RATE  CQ)’ .7X, ’THETA  DEG1./) 

STIME-STIME  +  CAPT 

IF  (  STINE  .GE.  DTI ME  )  STOP 

NUT-NUT  +1 

IF  (NUT  .EQ.  IPRNT)  GO  TO  3999 
GO  TO  77 
3999  TIME-STIME 

NUT  -  0 

CMATC 1.4) -X<1. 1) 

CMATC2. 4) -57. 29S78XXC  2. 1 )/( 1738. -X( 1 . 1 ) ) 
CMAT(3.4)-X(3. 1) 

CNAT(4.4)-X(4, 1) 

C  PUT  THE  APPROPIATE  UALUES  IN  DISKFILES 
URITEC5. 4001 )  TIME.  (CMATC  1 , 4) .  I -1 . 4) 

4001  FORMAT  C  TIME  -  * .F6.2. 1P4G15.5) 

URITEC1.200)  CMATC 1.4) 

URITEC2.200)  CMATC2.4) 

URITEC3.200)  CMATC3.4) 

WRITEC4.200)  CMATC4.4) 

GO  TO  22 
END 

C  MATMULT  IS  ROUTINE  TO  MULTIPLY  TWO  MATRIX  TOGETHER 
C 

SUBROUTINE  MATMUL  (AMATX, BMATX.CMAT.SIZE) 

DOUBLE  PRECISION  CMAT 
INTEGER  SIZE 

DIMENSION  AMATXC4.4) .BMATXC4.4) 

DIMENSION  CMAT  C4.4) 

DO  1  1-1. SIZE 
DO  2  J-l .SIZE 

2  CMAT ( I . J ) -0 . 0 

1  CONTINUE 

C  MULTIPLY  MATRICES 
DO  3  1-1. SIZE 
DO  4  K-l .SIZE 
DO  5  J-l. SIZE 
X-AMATX C I . J ) *BMATX C J . K ) 

5  CMAT ( I . K ) -CMAT ( I >  K ) +X 

4  CONTINUE 

3  CONTINUE 
RETURN 
END 

SUBROUTINE  INUSUB  (A.AI.RA.CA) 

DIMENSION  AC4.4) 

DIMENSION  IMC4. 4) . CMAT (4, 4 ) , AI (4,4) 

REAL  IM 
INTEGER  RA.CA 

C  INITIALIZE  ALL  MATRICIES 
DO  1  I  •  l.RA 

DO  2  J  •  l.RA 

AI ( I , J) -0. 0 


2 

1 

C 

c 


14 

c 

c 


imi.j)  . 

CONTINUE 


0.0 


FORM  THE  IDENTITY  MATRIX.  IM 
DO  14  I-l.RA 

inn.n-i.0 

CONTIfiJE 

initialize  the  inuerse  matrix 

START  BY  SETTING  AI  EQUAL  TO  A  MATRIX 
DO  3  I-l.RA 
DO  4  J-l.CA 
Ain.j)-A(i.j) 

CONTINUE 

IMJERT  THE  A  MATRIX 


FORWARD  ROW  REDUCTION  STARTS  HERE 
DO  S0  K-l.CA 
LINE  -  K 

9  FACTOR-A! (K.K) 

IF  (FACTOR  .EQ.  0.0)  (30  TO  2000 
LINE  -  K 
DO  71  I-l.RA 
AI (K. I )-AI (K. I ) /FACTOR 
IM(K.  I  )-IM(K.  D/FACTOR 
71  CONTINUE 

KPLUS1-K+1 
DO  20  I-KPLUS1  .CA 
FACTOR-AI ( I, K) 

DO  30  J-l.CA 

j  -AI  ( I .  J)-FACTORXAI  (K.  J) 

I ,  J)-FACTOR*IM(K.  J) 

CONTINUE 
CONTINUE 
CONTINUE 

REVERSE  ROW  REDUCTION  STARTS  HERE 
KSIZE-RA 
KSMl-KSIZE-1 
FACTOR-AI (KSM1.KSIZE) 

IF  (FACTOR  .EQ.  0.0)  GO  TO  330 
DO  310  J-l.RA 

J’ -AI  J)-FACTOR*AI  (KSIZE. 
J )- !  M  <  KSM1 .  J ) -FACTOR*  I M  ( KS !  ZE , 

KSMl-KSMl-1 

IF  (KSM1  .GT.  0)  GO  TO  90 
KSIZE-KSIZE  -1 
IF  (KSIZE  .GT.  1)  GO  TO  B9 
inuERSION  IS  COMPLETE  ASSIGN  AI  TO  IM 
DO  60  I-l.RA 
DO  61  J-l.CA 
AKI.J)-IM(I.J) 

CONTIMJE 
GO  TO  3002 

IF  (K  .EQ.  CA)  GO  TO  3000 
LINE-LINE+1 

.  IF  (LINE  .GT.  RA)  GO  TO  3000 
INTERCHANGE  LINES  OF  MATRIX 
LPl-LINE  +  1 
DO  320  J-l.CA 
CHAT  (K.J)-AICK.J) 

CMAT  (LPI.J)-IM(K.J) 

AKK.J)  -  AKLINE.J) 

IM(K.J)  •  INCLINE. J) 

AI(LINE,J)-CMAT  (K.J) 

IM(LINE.J)-CMAT  (LP1.J) 


30 

20 

S0 

C 

09 

90 


310 

330 


61 

60 

2000 


J) 

J) 
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Captain  Allan  T.  Reed  was  born  18  September  1953  in  Mountain 
Home,  Idaho.  He  graduated  in  1971  from  North  Miami  High  School  in 
Denver,  Indiana.  He  graduated  with  honors  from  Purdue  University  with 
a  Bachelor  of  Science  Degree  in  Aeronautical  and  Astronautical 
Engineering  in  1975.  Upon  graduation  he  received  a  commission  in  the 
USAF  as  an  ROTC  Distinguished  Graduate.  Captain  Reed  graduated  from 
pilot  training  from  Craig  AFB,  Alabama  in  1977  and  was  the  recipient 
of  the  ATC  Commander's  Trophy.  Initial  operational  experience  was  to 
Cannon  AFB,  New  Mexico  as  an  F-111D  pilot.  In  1980  he  was  selected  to 
attend  the  F-lll  Fighter  Weapons  School  and  upon  graduation  was 
selected  as  Distinguished  Graduate.  In  1982  Captain  Reed  served  as  an 
Instructor  Pilot  at  the  F-lll  Fighter  Weapons  School  at  Mountain  Home 
AFB,  Idaho.  In  1983  he  was  selected  for  the  Joint  Air  Force  Institute 
of  Techno logy/Test  Pilot  School  Program.  The  AFIT  course  work  was 
completed  in  June  1984.  Captain  Reed  graduated  from  the  USAF  Test 
Pilot  School  as  a  Distinguished  Graduate  in  June  1985. 


